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I.  Statement  of  Work 


During  the  1970’s  most  optimization  models  were  run  in  batch  mode  on  a  large 
mainframe  and  were  used  to  help  solve  some  type  of  planning  problem.  During  the 
1980’s  we  saw  more  demand  for  optimization  models  which  were  components  of  some 
real  time  system.  That  is,  decisions  were  being  made  sequentially  and  some 
optimization  model  was  communicating  with  an  on-line  database  and  was  being  used 
to  provide  information  for  real  time  decision  making.  In  the  1990’s,  clients  are 
demanding  real  time  systems  for  problems  which  were  formerly  solved  in  batch  mode. 

In  the  area  of  assigning  sailors  to  technical  schools  (both  C  and  A  Schools)  and 
assigning  sailors  to  job  billets,  the  Navy  is  already  developing  on-line  computer 
systems.  The  Job  Advertising  and  Selection  System  (JASS)  runs  on-line  and  allows  a 
career  counselor  on  a  ship  to  down  load  billet  information  from  a  server  located  in 
Washington  D.C.  A  version  of  JASS  which  contains  an  optimization  module  is  being 
developed  and  our  research  group  developed  the  optimization  software  which  is  being 
incorporated  into  the  system.  We  expect  additional  demand  for  high  speed 
optimization  software  in  the  area  of  Navy  personnel  assignment  and  our  work  on  the 
constrained  assignment  problem  is  related  to  this  demand. 

The  Navy  is  currently  developing  a  new  reservation  system  to  improve  the  assignment 
of  its  400,000  sailors  to  C  and  A  training  schools.  To  help  insure  the  success  of  the  new 
reservation  system,  modern  mathematical  models  and  computational  tools  are  needed 
to  help  generate  the  class  schedules  offered  by  the  various  training  schools.  The 
problems  are  members  of  the  class  NP-Hard  and  require  special  heuristic  algorithms 
for  real  time  application.  Our  work  on  efficient  algorithms  for  class  scheduling  is 
related  to  this  new  reservation  system. 

Other  Navy  applications  involve  real  time  targeting  of  cruise  missiles.  Generally  a 
missile  is  called  a  cruise  missile  if  its  speed  is  sub -sonic,  if  it  uses  a  built-in  global 
positioning  navigation  system  (GPS/INS),  and  if  its  range  is  at  least  several  hundred 
miles.  A  specific  mission  for  a  cruise  missile  is  programmed  by  specifying  a  sequence  of 
co-ordinates.  The  missile  uses  its  on— board  computational  facilities  and  its  GPS/INS 
system  to  guide  it  through  this  sequence  of  points.  At  present,  once  a  missile  is 
launched,  its  mission  can  not  be  modified.  However,  new  versions  are  under 
consideration  which  can  be  redirected  after  launch.  Finding  new  targets  involves 
on-line  optimization  and  our  work  on  efficient  primal  and  dual  algorithms  for 
generalized  networks  could  support  this  effort. 
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II.  Navy  Personnel  Assignment 


After  several  years  of  down  sizing,  the  Navy  currently  has  approximately  400,000  sailors 
each  of  whom  is  assigned  to  a  job  billet  for  a  tour  of  duty  lasting  from  two  to  five  years. 
Near  the  end  of  a  tour  (during  the  last  six  months)  the  sailor  is  placed  into  an 
assignment  pool  of  sailors  who  are  eligible  for  a  new  assignment.  Sailors  assigned  to  sea 
billets  are  usually  rotated  to  shore  billets  and  those  stationed  on  the  shore  are  usually 
assigned  to  a  sea  billet.  Some  200  detailers  within  the  Bureau  of  Navy  Personnel  are 
responsible  for  making  the  actual  assignments.  Every  one  that  has  studied  this  system 
agrees  that  Navy  personnel  assignment  is  a  very  complex  process  which  involves 
trade-offs  among  conflicting  policies.  It  is  also  complicated  by  the  fact  that  there  are 
many  more  billets  than  sailors. 

Over  the  last  decade,  operations  research  analysts  at  the  Navy  Personnel  Research  and 
Development  Center  have  developed  several  optimization  models  for  various  studies. 
Many  of  these  optimization  models  are  binary  integer  programming  problems  with 
special  structure.  The  underlying  structure  which  occurs  most  frequently  is  that  of  an 
assignment  problem  with  a  few  linear  side  constraints.  The  side  constraints  generally 
involve  a  budget  constraint  on  PCS  (permanent  change  of  station)  cost  and  may  involve 
quotas  on  seats  at  technical  schools. 

Our  research  team  has  developed  a  special  algorithm  for  the  constrained  assignment 
problem.  The  algorithm  is  based  on  branch— and— bound  and  uses  a  Lagrangean 
relaxation  to  improve  the  lower  bound.  A  complete  description  of  this  algorithm 
appears  in  Appendix  A. 


III.  Class  Scheduling  For  Navy  Training  Schools 


Before  assuming  a  new  assignment,  a  sailor  frequently  is  required  to  complete  one  or 
more  Navy  training  classes.  Experienced  sailors  (E4’s  through  E9’s)  attend  advanced 
skills  training  courses  while  new  recruits  (El’s,  E2’s,  and  E3’s)  take  a  set  of  basic 
courses.  The  advanced  courses  last  from  one  week  to  over  six  months  and  are  offered  by 
Navy  C  Schools.  The  basic  courses  are  generally  shorter  and  are  offered  by  Navy  A 
Schools.  Over  2000  different  courses  are  taught  each  year  at  a  cost  of  over  $1.3  billion. 

Navy  School  managers  develop  annual  class  offerings  which  are  then  used  by  the 
detailers  in  making  billet  assignments.  Until  recently,  the  class  schedules  were 
prepared  manually  by  personnel  at  each  School.  Due  to  budget  and  instructor 
reductions,  it  is  important  that  these  Navy  Schools  operate  as  efficiently  as  possible 
making  optimal  use  of  their  resources  (facilities  and  instructors). 

For  C— Schools  the  complicating  feature  is  that  many  sailors  who  attend  the  School 
need  several  different  courses.  A  sailor  only  attends  one  course  during  a  n  week  period, 
then  he/she  would  go  to  another  course.  Ideally  if  a  sailor  needed  four  courses,  then 
they  would  be  scheduled  back-to-back.  This  would  reduce  the  idle  time  between 
courses. 

The  A- School  courses  are  team  taught  with  different  numbers  of  instructors  required 
on  different  days  of  a  given  course.  The  objective  is  to  develop  a  feasible  schedule  that 
uses  the  smallest  number  of  instructors  on  the  busiest  day  of  the  year. 

Algorithms  for  these  scheduling  problems  have  been  developed  and  tested  on  real 
data.  The  results  of  our  study  may  be  found  in  Appendix  B. 


IV.  Generalized  Networks 


There  are  many  applications  in  the  area  of  network  flows  where  the  flow  either 
increases  or  decreases  as  it  passes  along  an  arc.  In  the  area  of  Navy  unit  personnel 
readiness,  promotions  and  attrition  can  be  modelled  using  this  feature. 

We  developed  a  new  optimizer  for  this  model  that  has  both  a  primal  and  dual  simplex 
implementation.  In  empirical  tests,  we  found  that  the  new  optimizer  is  approximately 
twenty  times  faster  than  CPLEX  3.0  on  moderate  sized  problems.  A  summary  of  our 
investigation  may  be  found  in  Appendix  C. 
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ABSTRACT 


This  manuscript  presents  a  truncated  branch -and— bound  algorithm  to  obtain  a 
near  optimal  solution  for  the  constrained  assignment  problem  in  which  there  are  only  a 
few  side  constraints.  At  each  node  of  the  branch-and— bound  tree  a  lower  bound  is 
obtained  by  solving  a  singly  constrained  assignment  problem.  If  needed,  Lagrangean 
relaxation  theory  is  applied  in  an  attempt  to  improve  this  lower  bound.  A  specialized 
branching  rule  is  developed  which  exploits  the  requirement  that  every  man  be  assigned 
to  some  job.  A  software  implementation  of  the  algorithm  has  been  tested  on  problems 
with  five  side  constraints  and  up  to  75,000  binary  variables.  Solutions  guaranteed  to  be 
within  10%  of  an  optimum  were  obtained  for  these  75,000  variable  problems  in  from 
two  to  twenty  minutes  of  CPU  time  on  a  Dec  Alpha  workstation.  The  behavior  of  the 
algorithm  for  various  problem  characteristics  is  also  studied.  This  includes  the  tightness 
of  the  side  constraints,  the  stopping  criteria,  and  the  effect  when  the  problems  are  un¬ 
balanced  having  more  jobs  than  men. 
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I.  INTRODUCTION 


The  constrained  assignment  problem  is  to  determine  a  least  cost  assignment  of  m 
people  to  n  jobs  such  that  an  additional  set  of  constraints  is  satisfied.  This  model  is  a 
binaiy  linear  program  and  may  be  stated  mathematically  as  follows: 


minimize  Cj-x- 

(1) 

(i.j)£A 

subject  to  ^  X”  =  1  ,  i  =  l,...,m 

(2) 

j:(i,j)eA 

X  xij  -  1  *  j =  1>-’n 

(3) 

i:(i,j)GA 

xij  e  {0,1},  all  (i,j)  e  A  (4) 

Z.  ^ijxij  -  *  k=l,...,s  (5) 

(ij)SA 


where  x-  =  1  implies  that  person  i  is  assigned  to  job  j  at  cost  of  c- ,  dH  denotes  the  coef¬ 
ficient  of  Xy  in  the  kth  side  constraint,  r k  denotes  the  right  -  hand  -  side  for  the  kth  side 
constraint,  and  A  is  the  set  corresponding  to  the  feasible  assignments.  Note  that  the 
problem  allows  for  more  jobs  than  people.  Many  practical  problems  have  this  feature. 
It  also  allows  for  the  case  in  which  the  number  of  people  exceeds  the  number  of  jobs. 
For  this  case,  one  simply  reverses  the  definition  of  people  and  jobs. 

The  problem  (1)  -(4)  with  m=n  is  the  classic  sparse  assignment  problem  which  is 
also  known  as  the  bipartite  matching  problem.  This  special  model  is  a  member  of  the 


class  P  and  there  are  excellent  algorithms  and  software  implementations  of  these  algo¬ 
rithms  available.  We  use  one  of  these  software  implementations  as  a  subroutine  in  this 
investigation. 

Since  (l)-(5)  is  a  binary  linear  program,  all  the  literature  on  integer  program¬ 
ming  applies  (see  Geoffrion  and  Marsten  [1972],  Salkin  [1974],  Parker  and  Rardin 
[1988],  Nemhauser  and  Wolsey  [1988]).  A  special  case  in  which  the  side  constraints 
have  the  generalized  upper  bound  (GUB)  structure  has  been  studied  by  Ali,  Kenning- 
ton,  and  Liang  [1993].  A  relaxation/decomposition  procedure  that  involves  solving  a 
series  of  pure  assignment  problems  is  used  successfully.  Ball,  Derigs,  Hilbrand,  and 
Metz  [1990]  also  present  an  algorithm  for  the  matching  problem  with  generalized  up¬ 
per  bound  side  constraints.  Another  special  case  for  s=l  and  m=n  has  been  studied  by 
Gupta  and  Sharma  [1981],  Aggarwal  [1985],  Bryson  [1991],  Kennington  and  Moham- 
madi  [1994, 1995].  The  only  specialized  algorithm  for  (l)-(5)  is  the  two  phase  proce¬ 
dure  of  Mazzola  and  Neebe  [1986].  The  first  phase  uses  subgradient  optimization  to 
obtain  an  advanced  startfor  the  branch-  and -bound  method  used  in  the  second  phase. 

This  research  was  motivated  by  our  work  with  the  Navy  Personnel  Research  and 
Development  Center  located  in  San  Diego.  This  general  model  often  appears  in  vari¬ 
ous  analytical  studies  which  involve  the  assignment  of  sailors  to  billets.  Most  of  these 
applications  are  for  the  unbalanced  model  in  which  there  are  more  jobs  (billets)  than 
sailors.  The  Navy  has  approximately  400,000  sailors  which  are  periodically  reassigned 
to  different  billets.  Sailors  rotate  between  sea  billets  and  shore  billets  with  various  re¬ 
strictions  on  these  assignments.  One  important  restriction  involves  the  cost  of  moving  a 
sailor  and  his/her  family  from  one  location  to  another.  This  is  referred  to  as  the  cost  for 
a  permanent  change  of  station  and  the  corresponding  side  constraint  is  known  as  the 


PCS  budget  constraint.  Fleet  balance  between  the  Atlantic  and  Pacific  Fleets  is  also 
important  and  leads  to  other  side  restrictions.  The  objective  of  this  study  was  to  develop 
a  truncated  branch- and-bound  algorithm  to  solve  the  constrained  assignment  prob¬ 
lem  and  to  provide  an  empirical  analysis  of  this  algorithm  on  a  variety  of  assignment 
problems  having  only  a  few  side  constraints. 


II.  THE  GENERIC  ALGORITHM 

Consider  the  problem  P(S)  s  min{  cx:  x  G  S}.  Using  the  terminology  of  Geoffrion 
and  Marsten  [1972],  v[P(S)]  denotes  the  optimal  objective  function  value  for  P(S), 
P  (S)  denotes  a  relaxation  of  P(S),  and  CL  denotes  the  candidate  list.  The  generic 
branch— and -bound  algorithm  used  in  this  investigation  may  be  stated  as  follows: 

Input: 

1.  The  problem,  P(S). 

Output: 

1.  The  solution  vector,  x* . 

2.  The  objective  value  corresponding  to  x' ,  v\(v‘  =  oo  implies  that  S  =  <!>.) 

Procedure  BAB; 

Begin 

initialize: 

CL:=  {P(S)},  v':=  oo; 
while  CL  *  <I>  do 

comment:  select  a  candidate  problem  for  analysis, 
select  P(U)  e  CL,  CL:=  CL\  {P(U)>; 
if  P(U)  has  a  feasible  solution,  then 
ifv[P(U)]  <  v’,  then 
let  x  be  an  optimum  for  P  (U); 
if  x  e  S,  then 
x'  :=  x  , v’:=  cx; 
else 


« 
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apply  a  heuristic  to  x  in  an  attempt  to  create  x  such  that  x  6  S  and 
cx  <  v'; 

if  successful,  then  x’:=  x,  v‘:=  cx; 
comment:  branching 

create  Ui,U2, Up such  that  U1UU2U  ...uUp  =  U  andUjOUj  =  O 
for  all  i*  j  e  {1, 2, p}; 

CL:=  CLu{P(Ui),  P(U2), P(Up)>; 

end  if 
end  if 
end  if 
end  while 
end. 
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III.  THE  SPECIALIZED  ALGORITHM 

The  specialized  techniques  developed  for  the  model  (1)  —  (5)  are  presented  in 
this  section.  This  includes  a  relaxation,  a  branching  rule,  and  fathoming  rules  based 
upon  the  underlying  assignment  structure.  These  are  combined  to  form  a  new  trun¬ 
cated  exponential  algorithm  for  the  constrained  assignment  problem. 

3.1  The  Relaxation 

Let  D  denote  the  matrix  corresponding  to  the  coefficients  in  (5),  x  denote  the  vec¬ 
tor  corresponding  to  the  binary  decision  variables,  c  denote  the  vector  of  costs,  and  r 
denote  the  vector  of  right-hand— side  values  for  the  side  constraints.  Then  the 
constrained  assignment  problem  may  be  denoted  as  P(S)  where  S  =  {x:  (2),  (3),  (4), 
(5)}.Letldenoteavectorof  l’sand  S  =  {x:  (2), (3), (4),lDx  <.ir}.ThenP(S)isavalid 
relaxation  for  P(S).  Application  of  the  algorithm  in  Kennington  and  Mohammadi 
[1994a]  to  solve  the  singly  constrained  assignment  problem  will  yield  a  lower  bound  for 
P(S),  the  optimal  Lagrangean  multiplier  corresponding  to  the  constraint  IDx  <Llr,  and 
x  e  S.Ifx  eS,  thencx  is  an  upper  bound  for  P(S). 

/V 

Let  S  =  {x:  (2),  (3),  (4)}.  Recall  that  a  Lagrangean  dual  for  P(S)  is  the  problem 

A 

max  {L(a)  :a  >.  0}  where  L(a)  =  min  {cx+  a(Dx— r):  xe  S  }.  We  use  the  optimal 
Lagrangean  multiplier  and  x  from  the  singly  constrained  algorithm  to  form  an  ad¬ 
vanced  starting  value  for  a .  Let  y  denote  the  solution  for  L  (a).  Then  w  =  Dy-r  is  used 
to  modify  a  for  successive  steps.  A  limited  number  of  these  steps  will  be  performed  in 
this  algorithm. 


3.2  The  Branching  Rule 

Consider  any  node  in  the  branch— and— bound  tree.  If  the  relaxation  P(  S)  has  no 
feasible  solution,  then  this  node  may  be  fathomed.  Otherwise,  an  assignment  will  be 
used  to  create  the  branches  as  illustrated  in  Figure  1. 

Figure  1  here 

For  a  given  node  in  the  branch -and— bound  tree  let  F'=  {(i,j):  x-  =  1}  and  F°= 
{(i,j):  Xjj  =  0}.  Let  U  =  {x  £  S:  x^  =  1  for  all  (i,j)  £  F'and  x-  =  0  for  all  (i,j)  e  F0}.  The 
relaxation  solved  at  each  node  in  the  branch-and-bound  tree  is  P(U)  which  is  a  singly 
constrained  assignment  problem  with  some  assignments  fixed.  Let  x  6  U,  T  =  {(i,j) : 
X-  =  l,(i,j)  $  F1},  andt=  |T|.  Consider  the  t+lsubsetsofU  (Ui,U2,...,U,+1)  created 
in  the  following  manner: 

Ui  =  {  x  €  U:  xiiji  =  0,  (i^)  e  T},Ti  =  ^{(inji)}; 

U2  =  {  x  £  U:  xiji  =  1,  Xjj^  =  0,  (i2j2)  £  Ti},  T2  =  T-|\{(i2j2)}; 

u3  =  {  x  £  U:  xiijt  =  1,  Xy2  =  1,  xijj3  =  0,  (i3,j3)  €  T2},  T3  =  T2\{(i3,j3)}; 

Ut  —  {  x  ^  U-  xiJi  —  1 , ...,  xit_jt  l-l,  Xj^  —  0,  (it*  jt )  £  Tt_,}; 

Ul+1=  {  x  £  U :  Xjj  =  1  for  all  (i,j)  £  T}  .  Note  that,  U  =  U1UU2U  ...  uut+1and 
UinUi=  <£>  for  all  i  ^  j. 

1  j  j 

Consider  a  node  in  the  branch-and-bound  tree  having  fi  edges  fixed  at  1.  Then 
branching  from  this  node  will  produce  t+ 1  =  n-fi  + 1  new  candidate  problems.  From 
Figure  1  it  can  be  seen  that  the  last  node  (i.e.,  Ut+)  need  not  be  created  since  it  was  ex¬ 
amined  at  the  parent  node.  Therefore,  each  branching  produces  t  candidate  problems. 
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3.3  The  Candidate  List 

For  our  implementation  of  the  algorithm,  a  singly  constrained  assignment  prob¬ 
lem  will  be  applied  to  P(Uj)  and  the  results  placed  in  the  candidate  list  (i.e.,  a  problem 
is  solved  before  it  is  placed  in  the  candidate  list).  The  motivation  for  placing  solved 
problem  in  the  candidate  list  is  that  the  solution  for  P(Ui)  can  be  easily  modified  to  ob¬ 
tain  an  advanced  starting  solution  for  P(Ui+).  Therefore  solving  the  sequence  of  prob¬ 
lems  P(Ut),  P(U2), ...,  P(Ut)  should  require  only  a  moderate  amount  of  computational 
effort.  Hence,  each  entry  in  the  CL  consists  of  the  five  tuple  ( F°,  F’  x  ,  (3 ,  u)  where  x  6 
U;,  |3  <.  v[P(Uj)],  and  u  the  optimal  Lagrangean  dual  for  the  singly  constrained  prob¬ 
lem. 

3.4  Fathoming  Rules 

At  any  node  p  of  the  branch-and-bound  tree  let  U,  F1,  and  F°be  as  defined  in 

Section  3.2.  Let  M  =  {1, 2, 3,...,  m}\  {i:  (i,j)  £  F1},  then  the  following  rules  may  be  used 
to  fathom  a  node.  If 

X  dij  +  Z  min(dij  :  (i’J)  e  A  \  F°)  >  rk 

(ij)€F»  i€M 

for  any  k,  then  node  p  can  be  fathomed.  That  is,  no  selection  of  the  free  variables  will 
satisfy  the  kth  side  constraint.  If 

Z  Cij  +  Z  min(Cij  :  (*’•>)  e  A  \  F°)  >  V* 

(i,j)£Fl  i£M 

then  node  p  can  be  fathomed.  That  is,  no  selection  of  the  free  variables  will  result  in  a 
solution  superior  to  the  incumbent. 

If  min  {IDx:  x  e  U  }  >  lr,  then  node  p  can  be  fathomed.  That  is,  no  selection  of  the 
free  variables  will  satisfy  all  side  constraints  simultaneously.  Let  |3  be  the  best  lower 


bound  obtained  for  node  p  and  s  be  the  termination  tolerance.  We  will  fathom  node  p  if 
v*  —  (3  <.  ep .  Using  this  rule  with  e  =  0.1  results  in  a  solution  from  the  procedure  guaran¬ 
teed  tobewithinlO%  ofthe  optimum  ande  =  0.01  produces  a  solution  within  1%  of  an 
optimum.  This  mle  speeds  convergence  at  the  expense  of  an  exact  solution. 

3.5  The  Algorithm 

In  this  section  the  information  presented  in  Sections  3. 1 —3.4  is  used  to  construct  a 
truncated  exponential  algorithm. 

Input: 

1.  The  cost  vector,  c. 

2.  The  feasible  region  S. 

3.  The  set  of  (man,  job)  pairs  corresponding  to  eligible  assignments,  A. 

4.  Termination  tolerance,  e. 

5.  The  maximum  execution  time,  tmax. 

6.  The  maximum  number  of  Lagrangean  relaxations  to  be  solved  at  each  node, 
limit. 

Output: 

1.  The  solution  vector,  x’ . 

2.  The  objective  value  corresponding  to  x‘ ,  v\  (v'  =  oo  implies  that  the  problem  is 
infeasible.) 

Procedure  ASSIGN+s; 

Begin 

initialize: 

comment:  Node  1  in  the  Branch- and— Bound  tree. 
v':=  oo,  F°:=  $,  Fh= 


A- 12 


ASSIGNP1(P(S),  x, P, u); 
if  P  =  —  oo  then  terminate; 
ifx  e  S  thenx'  :=  x  ,  v‘:=  cx; 
else  LAGRANGE(x  ,  P ,  u); 
comment:  Tolerance  test  for  fathoming, 
ifv'  —  p  <.  ep  then  terminate; 

CL:={(P,F>,x,p,u)}; 
while  CL  *  $  do 

SELECT  A  PROBLEMCF^F1,  x,  P,  u); 

BRANCH(CL,  F°,Fl,  x,  P,  u); 

end  while 
end. 

procedure  ASSIGNP1(P(U),  x,  p ,  u); 

Begin 

initialize: 

P  :=  -oo; 

apply  the  ASSIGN+1  algorithm  (Kennington  and  Mohammadi  [1994])  to  P(U); 
if  P(U)  has  a  feasible  solution  then 
let  x  be  the  best  feasible  solution  found  for  P(U); 
let  P  be  the  best  lower  bound  found  for  P(U); 

let  u  be  the  optimal  Lagrangean  dual  for  the  singly  constrained  problem; 
end  if 
end. 
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procedure  SELECT  A  PROBLEM  (F^F1,  x,  p,  u); 

Begin 

select  (F°,Fl,  x,  p,  u)  €  CL,  CL:=  CL  \  (F°,Fl,  x,  p,  u); 

end. 

procedure  BRANCH(CL,  F°,  F1,  x ,  p ,  u); 

Begin 

initialize: 

t:=  0,  G:=  {(i,j):  xAj  =  1,  (i,j)  $  F1}; 

M:=  {1, 2, 3, m}\  {i:  (i,j)  e  F1}; 
while  G  *>*  <E>  do 

comment:  Fix  a  variable  at  zero. 

let  (itj 0  e  G,  G:=  G\  {(i^)},  F“:=  P  u  {(i^)}; 

ASSIGNP1(P(U),  x,  p ,  u); 
if  P  s*  — oo  then 

if  x  e  S  and  cx  <  v*  then  x'  :=  x  ,  v‘:=  cx; 
else  LAGRANGE(x  ,  p ,  u); 
comment:  Tolerance  test  for  fathoming, 
if  v*  -  p  >  ep  then  CL:=  CL  u  {(F0,F‘,  x,  p,  u)}; 
comment:  Permanently  assign  man  U  to  job  ji. 

M:=M\{U},  F':=  Pu 

F°:=  F°  u  {(i-i,j)  :  (ii,j)  e  Afor  all  j  }  u  { (i,ji) :  (i,ji)  s  Afor  all  i  }\  {(ii,ji)}; 
comment:  Assignment  polytope  feasibility  tests. 

Cjj  +  min(Cij :  (i,  j)  G  A)  >  v*  then  return; 

(i,j)€F>  ieM 


if  dji  +  min(dij :  OJ)  e  A)  >  rk  then  return 
(ij)eF'  iGM 

end  for 

end  if 

end  while 

end. 

Procedure  LAGRANGE(x  ,  |3 ,  u); 

Begin 

initialize: 

a:=  ui,  y:=x  ,t:=  1; 
while  (v* — (3  >  e|3  and  t  <.  limit)  do 
w:=  Dy-r; 
for  i=l,...,s 

if  wk  >  0,  then  ak:=  1.25afc; 
if  wt  <  0,  then  a*:  =  0.75a,.; 

end  for 

let  y  be  a  solution  for  L(a)  and  p  :=  max{p ,  v[L(a)]}; 
if  y  G  S  and  cy<  v*  then  x’  :=  y,  v‘:=  cy; 
t:=  t+1; 

end  while 


This  algorithm  exploits  the  structure  of  the  model  (l)-(5)  in  several  ways.  Per¬ 
manent  assignment  of  a  man  to  a  job  implies  that  all  other  variables  involving  this  man 
and  job  may  be  fixed  at  zero.  The  relaxation  is  a  singly  constrained  assignment  problem 
for  which  near  optimal  integer  solutions  can  be  obtained  using  the  results  in  Kenning- 
ton  and  Mohammadi  [1994].  Special  fathoming  rules  which  are  based  upon  the  assign¬ 
ment  polytope  are  used.  • 


IV.  EMPIRICAL  ANALYSIS 


The  specialized  algorithm  has  been  implemented  in  software  (called  AS¬ 
SIGN +s)  and  empirically  analyzed  on  an  Alpha  workstation  by  Digital  Equipment 
Corporation.  The  code  is  written  in  Fortran  and  uses  ASSIGN  + 1  (see  Kennington  and 
Mohammadi  [1994])  to  solve  the  singly  constrained  assignment  problems.  ASSIGN + 1 
is  an  implementation  of  the  Lagrangean  relaxation  algorithm  for  sparse  singly 
constrained  assignment  problems. 

A  test  problem  generator  was  developed  which  has  the  following  inputs:  (i)  the 

number  of  men,  (ii)  number  of  jobs  for  each  man,  (iii)  the  maximum  cost,  c,  (iv)  the 
number  of  side  constraints,  s,  and  (v)  the  side  constraint  multiplier,  z.  Both  the  costs 

and  the  side  constraint  coefficients  are  uniformly  distributed  over  the  range  (0,  c).  We 
randomly  generate  a  feasible  assignment,  x,  and  set  the  right— hand -side  of  the  side 

constraints,  r,  to  zDx.  Obviously,  as  z  becomes  smaller,  the  feasible  region  becomes 
smaller  and  for  sufficiently  small  z  {x:  (2),  (3),  (4),  and  (5)}  is  usually  empty. 

The  generator  was  used  to  generate  two  sets  of  400x400  problems  described  in 
Table  1.  As  Table  1  indicates,  the  problems  generally  become  more  difficult  as  z  be¬ 
comes  smaller  and  for  z  small  enough  the  feasible  region  is  empty.  For  all  runs,  the  stop¬ 
ping  criteria  used  is  £=10%  and  the  %  deviation  reported  in  column  8  gives  a  guarantee 
on  the  deviation  from  optimality.  All  times  are  CPU  time  and  exclude  the  time  for  both 
input  and  output.  The  run  with  problem  number  2  having  z=0.4  was  terminated  after 
the  candidate  list  grew  to  25,000  entries.  As  we  expected,  there  exists  problems  which 
cannot  be  solved  in  a  reasonable  amount  of  time  and  storage  using  this  approach.  Tight- 
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ly  constrained  problems  having  48,000  binary  variables  definitely  stretches  the  capabil¬ 
ity  of  this  software  implementation. 

Table  1  here 

Tables  2  and  3  give  our  empirical  results  with  30  randomly  generated  assignment 
problems  with  various  sizes  all  having  five  side  constraints.  For  all  of  the  problems 
tested  we  were  able  to  find  a  solution  guaranteed  to  be  within  10%  of  an  optimal  solu¬ 
tion.  The  six  smallest  problems  have  3,000  binary  variables.  Five  of  these  were  solved  in 
less  than  two  minutes  each  and  one  required  about  six  and  one -half  minutes.  The  six 
largest  problems  had  75,000  binary  variables  and  were  all  solved  in  less  than  twenty  one 
minutes  each.  The  most  difficult  problems  (300x300)  have  27,000  binary  variables.  Two 
of  these  six  problems  required  eighty  minutes  to  solve.  These  two  difficult  problems 
also  had  very  tight  side  constraints. 

Tables  2  and  3  here 

This  work  was  motivated  by  models  for  assigning  sailors  to  ships  and  for  this  ap¬ 
plication  the  number  of  jobs  always  exceeds  the  number  of  sailors  available.  Frequently 
the  j ob  list  covers  a  longer  period  than  the  list  of  available  sailors  which  produces  a  large 
imbalance  in  n  and  m.  Tables  4  and  5  present  our  empirical  results  from  solving  18  un¬ 
balanced  assignment  problems  with  five  side  constraints.  For  the  300x600  problem  with 
z=0.6  presented  in  Table  5,  the  run  was  terminated  due  to  candidate  size  limit.  For  all 
other  test  problems  we  were  able  to  obtain  a  solution  within  10%  of  an  optimum. 

Tables  4  and  5  here 

For  all  test  problems  we  search  for  a  solution  within  10%  of  an  optima.  To  study 
the  effect  of  the  tolerance  value  on  the  performance  of  the  algorithm  we  solved  two 
200x200  and  two  200x400  problems  with  different  tolerance  values.  Figure  2  indicates 


that,  as  expected,  a  decrease  in  the  tolerance  value  leads  to  an  increase  in  the  execution 
time.  For  all  four  problems,  a  point  was  reached  in  which  a  slight  decrease  in  the  toler¬ 
ance  resulted  in  a  large  increase  in  the  solution  time. 

Figure  2  here 


V.  SUMMARY  AND  CONCLUSIONS 

We  have  presented  a  truncated  exponential  algorithm  for  the  constrained  assign¬ 
ment  problem.  The  algorithm  is  applicable  for  both  balanced  and  unbalanced  assign¬ 
ment  problems  having  inequality  side  constraints.  The  algorithm  uses  a  specialized 
branching  rule  that  exploits  the  underlying  structure  of  the  problem.  Bounds  are  ob¬ 
tained  by  solving  a  singly  constrained  assignment  problem  followed  by  a  few  iterations 
with  a  Lagrangean  relaxation. 

We  present  empirical  results  for  both  balanced  and  unbalanced  problems  having 
five  side  constraints.  For  problems  having  75,000  binary  variables,  solutions  guaran¬ 
teed  to  be  within  10%  of  an  optima  were  obtained  in  less  than  twenty  one  minutes  on  a 
Dec  Alpha  workstation.  Our  analysis  indicated  that  as  the  side  constraints  become 
tighter  the  execution  time  and  number  of  branch— and— bound  nodes  increases.  For 
one  of  the  300x300  problems  having  27,000  arcs,  the  execution  time  increased  from 
about  one  minute  to  about  eighty  minutes  as  a  result  of  side  constraint  tightening.  Our 
analysis  also  indicates  that  the  performance  of  the  algorithm  on  unbalanced  problems 
is  generally  better  than  its  performance  for  the  balanced  problems  with  the  same  num¬ 
ber  of  binary  variables.  The  Navy  personnel  assignment  problems  which  motivated  this 
study  are  all  unbalanced  models. 

For  problems  of  this  type,  having  only  a  few  side  constraints,  we  believe  that  this  is 
the  current  best  algorithm  and  software  implementation  available.  Solutions  guaran¬ 
teed  to  be  within  10%  of  an  optimum  should  be  obtained  for  most  problems  having  few¬ 
er  than  five  side  constraints  and  fewer  than  20,000  arcs. 
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a.  Problem  #1  (200x200) 
with  1200  arcs 


b.  Problem  #2  (200x200) 
with  1200  arcs 


Time  (min.)  Time  (min.) 


c.  Problem  #3  (200x400)  d.  Problem  #4  (200x400) 

with  1200  arcs  with  1200  arcs 
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Table  1.  Empirical  results  from  the  algorithm  for  400x400  assignment  problems  with  five  side 

constraints  (48,000  columns  and  805  rows). 


Prob. 

z 

#  BAB 
Nodes 

#AP’s 

Solved 

LB 

UB 

Per  Cent 
Deviation 

Node  #  of 
Incumbent 

1.0 

189 

216 

0.71 

5.270 

5,309 

0.74 

13 

0.9 

244 

1,958 

4.74 

5,505 

5,768 

4.78 

22 

0.8 

274 

2,119 

3.96 

6,999 

7,271 

3.89 

78 

1 

0.7 

952 

8,288 

12.66 

10,991 

11,725 

6.68 

859 

0.6 

14.142 

mi 

166.00 

20.820 

22.713 

9.09 

13.943 

0.5 

4,214 

35,362 

45.39 

47,446 

50,343 

6.11 

4,157 

0.4 

1 

2 

0.00 

problem  has  no  feasible  solution 

1.0 

224 

503 

1.66 

5,370 

5,559 

3.52 

2 

0.9 

253 

1,801 

4.57 

5,830 

6,132 

5.23 

34 

0.8 

700 

5,434 

9.26 

7,815 

8,373 

7.13 

492 

2 

0.7 

3,114 

25,786 

42.93 

12,510 

13,649 

9.10 

2,938 

0.6 

3,094 

28,694 

40.92 

23,606 

25,188 

6.70 

2,871 

0.5 

3,505 

29,642 

31.45 

50,939 

54,690 

736 

3,340 

25,0001 

79,160 

53.37 

148,825 

no  feasible  solution  obtained 

0.3 

1 

2 

0.0 

problem  has  no  feasible  solution 

i  terminated  due  to  candidate  list  size  limit. 


Table  2.  Empirical  results  with  problem  set  1  using  the  algorithm. 
_ (All  problems  have  five  side  constraints.) _ 


Table  3.  Empirical  results  with  problem  set  2  using  the  algorithm. 
_ (All  problems  have  five  side  constraints.) _ 
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Table  4.  Empirical  results  with  problem  set  3  using  the  algorithm. 

(All  problems  have  live  side  constraints.) _ 
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Table  5.  Empirical  results  with  problem  set  4  using  the  algorithm. 

(All  problems  have  five  side  constraints.) _ 
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Abstract 


The  problem  of  developing  good  schedules  for  Navy  C— Schools  has  been  modeled  as  a 
combinatorial  optimization  problem.  The  only  complicating  feature  of  the  problem  is  that 
classes  must  be  grouped  together  into  sequences  known  as  pipelines.  An  ideal  schedule  will 
have  all  classes  in  a  pipeline  scheduled  in  consecutive  weeks.  The  objective  is  to  eliminate  the 
nonproductive  time  spent  by  sailors  at  C— Schools  who  are  waiting  for  the  next  class  in  a  pipe¬ 
line.  In  this  investigation,  five  algorithms  were  specialized  for  this  problem  and  empirically 
analyzed  on  a  set  of  actual  scheduling  problems.  The  algorithms  evaluated  include  simulated 
annealing,  a  greedy  heuristic,  tabu  search,  a  genetic  algorithm,  and  evolutionary  program¬ 
ming.  In  empirical  testing  the  greedy  heuristic  was  found  to  be  best  for  this  application. 

In  addition,  an  implicit  enumeration  procedure  for  this  problem  was  also  developed. 
The  implicit  enumeration  algorithm  uses  the  output  from  the  greedy  algorithm  as  the  initial 
upper  bound  and  develops  a  lower  bound  based  on  the  problem  input  data.  If  the  bounds  are 
equal,  then  an  optimal  schedule  has  been  found;  otherwise,  an  implicit  enumeration  proce¬ 
dure  is  initiated  which  if  allowed  to  run  long  enough  will  eventually  find  an  optimal  schedule. 
On  five  of  the  six  test  problems,  a  provable  optimum  was  obtained.  The  sixth  problem  remains 
unsolved. 

The  best  ideas  for  Navy  C- School  scheduling  were  extended  for  basic  training  courses 
which  are  offered  in  Navy  A-Schools.  In  empirical  tests,  we  found  that  our  heuristic  algorithm 
worked  very  well  on  the  A-School  test  problems. 

Acknowledgment 

This  work  was  supported  by  the  Navy  Personnel  Research  &  Development  Center  under  the  auspices  of  the  U.  S. 
Army  Research  Office  Scientific  Services  Program  administered  by  Battelle  ( Delivery  Order  1110,  Contract 
Number  DAAL03-91-C-0034)  and  by  the  Office  of  Naval  Research  under  Contract  Number  N000 14-95- 1-0645. 


I.  INTRODUCTION 


The  Navy  currently  has  approximately  400,000  sailors  each  of  whom  is  assigned  to  a  job 
billet  for  a  specific  tour  of  duty  lasting  from  two  to  five  years.  After  the  tour  is  completed,  a 
sailor  is  assigned  to  a  new  billet  for  another  tour  of  duty.  Sailors  assigned  to  sea  billets  are 
usually  rotated  to  shore  billets  and  vice  versa.  Each  month,  the  Navy  assigns  thousands  of  new 
recruits  and  rotating  sailors  to  vacant  billets  (jobs).  The  assignments  are  made  by  some  200 
detailers  within  the  Personnel  Assignment  Division  of  the  Bureau  of  Navy  Personnel.  Enlisted 
personnel  assignment  is  a  complex  process  which  often  involves  trade-offs  among  several 
conflicting  policies.  An  excellent  description  of  the  Navy5  s  system  may  be  found  in  Blanco  and 
Hillery  [1994]. 

Before  assuming  a  new  assignment,  a  sailor  frequently  is  required  to  complete  one  or 
more  Navy  training  classes.  Experienced  personnel  attend  advanced  skills  training  courses 
while  new  recruits  take  a  set  of  basic  courses.  The  advanced  courses  are  offered  by  Navy  C— 
Schools  and  the  basic  courses  are  offered  by  Navy  A— Schools.  Multiple  offerings  of  over  2000 
different  courses  must  be  taught  each  year  with  a  course  lasting  from  one  week  to  over  six 
months.  The  Navy’s  annual  budget  for  training  exceeds  $1.3  billion. 

Based  on  the  forecast  demand  for  courses,  the  C— School  and  A— School  managers  de¬ 
velop  a  schedule  of  class  offerings  for  use  by  the  detailers.  Until  recently  these  schedules  were 
prepared  manually  by  personnel  at  each  of  the  Schools.  The  creation  of  an  optimal  schedule 
for  either  a  C- School  or  an  A- School  is  a  difficult  combinatorial  optimization  problem.  The 
objective  of  this  study  is  to  develop  and  empirically  evaluate  algorithms  for  creating  good  class 
schedules  for  both  types  of  Navy  Schools. 


II.  NAVY  C- SCHOOLS 


The  assignment  to  a  new  billet  may  require  that  the  sailor  update  his/her  skills  which  is 
accomplished  by  taking  courses  in  the  Skill  Progression  Training  Program.  These  courses  are 
taught  at  several  campuses  in  the  continental  U.S.A.  with  the  largest  campuses  being  at  the 
Fleet  Training  Center  in  Norfolk  and  at  the  Service  School  Command  Headquarters  in  San 
Diego.  The  organizations  which  offer  these  courses  are  called  C-  Schools.  During  FY95  these 
C— Schools  offered  several  thousand  different  courses  to  approximately  60,000  enlisted  per¬ 
sonnel. 

As  the  military  draw  down  proceeds  and  budgets  for  skills  training  is  decreased,  it  is 
imperative  that  these  Schools  operate  as  efficiently  as  possible.  An  important  problem  faced 
by  the  management  of  these  schools  is  the  development  of  a  good  class  schedule.  In  this  section 
we  describe  the  pipeline  scheduling  problem,  present  a  set  of  algorithms  for  this  problem, 
provide  an  empirical  evaluation  of  our  algorithms,  and  finally  give  our  recommendation  for 
solving  this  problem. 


2.1  The  Pipeline  Scheduling  Problem 

The  Service  School  Command  at  the  San  Diego  Naval  Station  offers  four  courses  in  the 
area  of  maintenance  of  electronic  communication  equipment  for  the  FFG-7  ship  class.  A 
sailor  who  successfully  completes  all  four  courses  is  called  a  Small  Combatant  Communications 
Subsystem  Technician.  The  courses  are  as  follows: 

(i)  HFSYS  -  (High  Frequency  Systems  -  12  weeks), 

(ii)  WSC3  —  (AN/WSC-3  Communications  Sets  Maintenance  -  12  weeks), 

(iii)  NAVMACS  -  (Organization  Level  Maintenance  on  the  Navy  Modular  Automated  Com- 


munications  System  -  12  weeks),  and 

(iv)  VRC46  —  (AN/VRC— 46  Radio/Transmitter  Set  Maintenance  —  2  weeks). 

During  1995  this  School  offered  HFSYS  and  VRC46  eight  times,  WSC3  nine  times,  and 
NAVMACS  sixteen  times.  A  class,  as  distinguished  from  a  course,  is  a  specific  offering  (con¬ 
vening  date)  of  a  course.  Hence,  during  the  one  year  planning  horizon,  there  will  be  eight 
classes  of  courses  HFSYS  and  VRC46,  nine  classes  of  WSC3,  and  sixteen  classes  of  NAV¬ 
MACS.  Since  a  sailor  can  only  take  one  class  at  a  time,  a  sailor  who  needs  all  four  courses  will 
be  assigned  to  this  campus  for  a  minimum  of  thirty  eight  weeks. 

This  School  runs  two  shifts  (a  day  shift  and  an  evening  shift)  for  courses  HFSYS,  WSC3, 
and  NAVMACS,  and  a  day  shift  only  for  VRC46.  Many  sailors  come  for  all  four  courses  and 
proceed  through  four  classes  grouped  together  into  what  is  known  as  a  pipeline  or  pipe.  That 
is,  a  pipe  is  a  sequence  of  classes  that  are  scheduled  so  that  a  sailor  can  take  all  four  courses 
before  departing  for  the  next  assignment.  The  four  courses  in  a  pipe  can  be  taken  in  any  order 
and  ideally  they  will  be  offered  back-to-back. 

Classes  are  generally  team  taught  and  there  may  be  several  teams  available  to  teach  the 
same  course.  When  there  are  multiple  teams,  they  teach  during  different  shifts  (either  the  day 
shift,  the  evening  shift,  or  the  night  shift).  Hence,  WSC3  shift  1  and  WSC3  shift  2  refer  to  the 
same  course  (WSC3)  being  taught  by  different  teams  at  different  times  of  the  day.  For  this 
study  the  terms  team  and  shift  will  be  used  interchangeably. 

Since  many  courses  have  multiple  instructors,  it  is  sometimes  possible  for  a  given  team 
to  begin  a  new  class  before  a  previous  class  is  finished.  The  minimum  time  required  between 
the  starting  times  of  two  classes  is  called  the  offset.  For  example,  suppose  course  A  lasting  10 
weeks  has  a  two  week  offset.  Then  classes  could  begin  in  weeks  one,  three,  five,  seven,  et  cet¬ 
era.  This  partial  schedule  is  illustrated  in  Figure  1. 


Figure  1.  About  Here 

Ideally  the  classes  in  a  pipe  will  be  offered  back— to— back  so  that  a  sailor  has  no  idle 
weeks  at  the  C-  School.  A  period  of  one  or  more  weeks  in  a  pipe  in  which  no  class  is  offered  is 
called  a  gap  in  the  pipe.  Pipe  1  illustrated  in  Figure  2  has  no  gaps  while  pipe  2  has  gaps  totaling 
four  weeks.  The  pipes  in  a  perfect  schedule  will  have  no  gaps. 


Figure  2.  About  Here 


These  ideas  are  further  illustrated  by  examining  the  problem  data  for  the  FFG— 7  prob¬ 
lem  presented  in  Table  1.  Line  1  indicates  that  this  schedule  is  for  the  fiscal  year  1995  (1  Octo¬ 
ber  1994  -  30  September  1995),  and  is  a  one  year  schedule.  A  one  year  schedule  corresponds 
to  50  weeks  with  a  two  week  break  for  Christmas.  A  two  year  schedule  corresponds  to  100 
weeks  and  all  schedules  are  for  at  most  two  years. 


Table  1.  About  Here 


Line  3  indicates  that  the  schedule  is  composed  of  eight  pipes  each  of  which  has  four 
courses.  The  names  of  the  four  courses  are  given  on  lines  4, 8, 12,  and  16.  Course  A  lasts  12 
weeks  and  there  are  two  teams  (shifts)  available  to  teach  it.  The  first  shift  will  offer  this  course 
four  times  during  the  year  with  an  offset  equal  to  the  course  length.  That  is,  the  second  offering 
of  course  A  using  shift  1  cannot  begin  until  the  completion  of  the  first  offering.  For  course  B, 
shift  1  will  offer  five  classes  during  the  year  while  shift  2  will  only  offer  four.  This  gives  a  total  of 
nine  class  offerings  for  the  eight  pipes.  Hence,  one  of  the  classes  will  not  appear  in  any  pipe. 
An  offset  of  eight  for  course  B  shift  1  indicates  that  if  a  class  begins  in  week  t  another  class  could 
begin  in  week  t+8.  Hence,  in  weeks  t+8,  t+9,  t+10,  and  t+ 11  there  wouldbe  two  classes  being 


taught  by  shift  1 .  For  course  C  there  will  be  a  total  of  16  classes  of  which  only  eight  will  appear 
in  pipes. 

Based  upon  its  characteristics,  a  schedule  is  assigned  a  nonnegative  score.  Penalty  points 
are  assessed  for  undesirable  characteristics  with  a  perfect  schedule  having  a  score  of  zero. 
Gaps  in  pipes  result  in  penalty  points  equal  to  the  total  number  of  weeks  in  the  gaps.  A  class 
which  extends  beyond  the  planning  period  is  assessed  penalty  points  equal  to  the  number  of 
weeks  that  it  extends  beyond  the  planning  horizon.  For  example  a  class  that  was  completed  at 
the  end  of  week  56  when  the  planning  horizon  is  50  weeks  would  be  assessed  a  penalty  of  six 
points.  Each  pipe  is  assigned  a  minimum  target  week  for  its  start  date.  A  pipe  is  assessed  one 
penalty  point  for  each  week  of  time  that  it  begins  prior  to  its  minimum  target  week.  Pipes  may 

begin  on  or  after  the  target  week  with  no  penalty. 

The  minimum  pipeline  target  for  pipe  p  (Ep)  is  based  on  the  number  of  weeks  in  the  plan¬ 
ning  horizon  (W)  and  the  number  of  pipes  (P)  in  the  schedule.  This  relationship  is  Ep  = 
L(p-l)W/(2P)J+  1.  ForFFG-7,  the  pipeline  targets  for  each  ofthe  eight  pipelines  are  1,4, 7, 
10, 13, 16, 19,  and  22.  A  feasible  schedule  having  a  score  of  19  appears  in  Figure  3.  The  eight 
pipes  are  denoted  by  the  lower  case  letters  a,  b,  ...,h.  Pipe  a  consists  of  the  classes  Ala,  Die, 
Cld,  Bid.  The  0’s  in  Figure  3  refer  to  classes  which  do  not  appear  in  any  pipe.  The  detailed 
calculations  of  this  score  may  be  found  in  Tables  2  and  3.  The  six,  three,  and  ten  in  Table  2  are 
due  to  the  three  classes  which  extend  beyond  the  50  week  planning  horizon. 


Figure  3,  Table  2,  and  Table  3.  About  Here 


2.2  Heuristic  Algorithms  for  the  Pipeline  Scheduling  Problem 

Heuristic  procedures  for  discrete  optimization  problems  are  usually  based  upon  one  of 
three  ideas:  (i)  greedy  heuristics,  (ii)  interchange  heuristics,  or  (iii)  truncated  exponential  pro¬ 
cedures.  This  section  presents,  a  greedy  heuristic  and  several  interchange  heuristics  which 


have  been  specialized  for  the  pipeline  scheduling  problem.  The  interchange  heuristics  include 
simulated  annealing,  tabu  search,  a  genetic  algorithm,  and  evolutionary  programming.  Some 
general  information  about  all  of  these  techniques  can  be  found  in  Parker  and  Rardin  [1988], 
and  Nemhauser  and  Wolsey  [1988, 1989]. 

2.2.1  Simulated  Annealing 

The  term  simulated  annealing  comes  from  the  thermodynamics  analogy  with  the  way 
that  metals  cool  or  anneal  (see  Press,  Flannery,  Teukolsky,  and  Vetterling  [1989]).  If  a  liquid 
metal  is  cooled  slowly,  then  the  atoms  form  a  crystal  which  has  a  minimum  energy  state.  If  a 
metal  is  cooled  too  quickly  or  quenched,  then  it  may  end  up  in  a  polycrystalline  state  having  a 
higher  energy  level.  The  idea  behind  simulated  annealing  is  to  use  slow  cooling  in  an  attempt 
to  obtain  a  low  energy  state  (objective  value). 

A  software  implementation  of  the  simulated  annealing  algorithm  specialized  for  the 
pipeline  scheduling  problem  is  called  Pipeline  Assistant.  It  begins  with  a  set  of  equally  spaced 
start  times  for  each  class  and  a  random  assignment  of  classes  to  pipes.  This  starting  schedule 
may  violate  both  the  offset  restrictions  and  the  restriction  that  a  sailor  cannot  be  in  two  classes 
simultaneously.  A  very  large  penalty  is  assessed  for  these  constraint  violations.  At  each  itera¬ 
tion  a  neighboring  schedule  is  generated  at  random.  Generating  a  neighbor  will  involve  either 
moving  a  single  class  one  week  earlier  or  one  week  later,  or  switching  a  pair  of  classes  of  the 
same  course  between  two  pipes. 

If  the  neighbor  has  a  lower  score  than  the  current  schedule,  then  the  neighbor  becomes 
the  current  schedule  and  the  process  is  repeated.  If  the  neighbor’s  score  is  worse  than  that  of 
the  current  schedule,  then  a  random  number  is  used  to  determine  if  the  neighbor  becomes  the 
current  schedule.  Suppose  the  current  schedule  has  a  score  of  Sc  and  the  neighbor  has  a  score 
of  Sn  with  Sc  <  Sn.  Let  R  be  a  random  number  from  a  Uniform  [0,1]  distribution,  andletTbea 
parameter  known  as  the  temperature.  If  R  <  [1  —  (Sn  —  Sc)  /T],  then  the  neighbor  becomes 
the  new  current  schedule.  Initially  T  is  set  to  5%  of  the  score  of  the  initial  schedule.  As  the 
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iterations  proceed,  the  temperature  is  systematically  lowered  by  replacing  T  with  0.8T.  The 
procedure  searches  for  some  finite  number  of  iterations  and  retains  the  best  schedule  found 
during  the  search.  In  Pipeline  Assistant  this  is  repeated  three  times  and  these  three  schedules 
are  reported  to  the  user.  There  is  no  guarantee  that  any  of  these  schedules  will  be  feasible. 

2.2.2  A  Greedy  Heuristic  With  Local  Improvements 

The  greedy  heuristic  develops  a  schedule  similar  to  the  way  a  human  would  construct  a 
schedule  manually.  Initially  the  courses  are  given  a  priority  ranking.  For  a  four  course  pipe¬ 
line,  course  A  is  given  priority  1,  course  B  is  given  priority  2,  and  so  forth.  The  pipes  all  begin  at 
their  minimum  target  times  and  are  constructed  one  class  at  a  time.  At  each  iteration,  the 
shortest  incomplete  pipe  is  selected  for  a  class  assignment.  The  first  available  class  associated 
with  the  unassigned  course  having  highest  priority  is  assigned  to  this  pipe  at  the  earliest  pos¬ 
sible  time.  This  strategy  attempts  to  construct  pipes  having  no  gaps.  After  all  pipes  have  been 
constructed,  the  other  classes  are  scheduled  one  at  a  time  at  the  earliest  time  possible.  At  this 
stage  a  feasible  schedule  has  been  developed. 

Each  pipe  having  a  gap  of  at  least  one  week  is  examined  in  an  attempt  to  start  the  pipe 
later  and  thus  reduce  the  gap.  Finally,  for  those  pipes  in  which  the  last  class  begins  after  the 
planning  horizon,  an  attempt  is  made  to  move  this  class  to  the  first  position  in  the  pipe.  Fewer 
penalty  points  will  be  incurred  for  beginning  this  pipe  prior  to  its  target  week  than  are  currently 
being  assessed  for  beginning  after  the  planning  horizon.  This  results  in  one  candidate  sched¬ 
ule. 

Other  candidate  schedules  are  obtained  by  changing  the  course  priorities.  The  greedy 
algorithm  tries  a  large  number  of  possible  priorities  and  saves  the  best  schedule  found.  For  a 
problem  having  n  courses,  min(10K,  n!)  possible  schedules  will  be  developed.  Some  theoreti¬ 
cal  results  concerning  the  optimality  of  greedy  algorithms  on  ordering  problems  may  be  found 
in  Dechter  and  Dechter  [1989]. 
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2.2.3  The  Tabu  Search  Algorithm 

Tabu  search  is  a  heuristic  technique  that  is  based  upon  selected  concepts  from  the  field  of 
artificial  intelligence.  It  is  a  flexible  method  for  guiding  a  search  over  the  space  of  solutions  in 
an  attempt  to  discover  a  high  quality  solution.  Tabu  search  uses  three  types  of  rules  that  at¬ 
tempt  to  prevent  the  search  from  being  trapped  at  an  inferior  local  optimum. 

The  first  set  of  rules  guides  the  local  search.  Let  w  be  a  given  number  of  weeks  that  a  class 
is  allowed  to  move,  either  forward  or  backward.  For  a  given  schedule  x  and  a  given  class  c,  the 
neighborhood  N(x,c),  is  defined  as  the  set  of  all  schedules  formed  from  xby  moving  c  at  most  w 
weeks.  A  neighborhood  N(x,c)  is  said  to  be  tabu  if  class  c  has  been  moved  during  the  most 
recent  q  iterations.  The  systematic  use  of  this  rule  provides  the  tabu  search  with  a  short  term 
memory  that  prevents  the  local  search  from  revisiting  the  same  neighborhood.  The  “most  im¬ 
proving”  selection  rule  would  select  the  best  move  over  all  N(x,c)  for  all  classes  c.  For  this  ap¬ 
plication  the  most  improving  rule  is  numerically  expensive  and  requires  the  evaluation  of  each 
move  generated  by  each  (course,  shift,  class,  week)  combination.  In  our  tabu  search  proce¬ 
dure,  we  use  a  hybrid  of  the  most  improving  rule  and  a  Monte  Carlo  method.  Specifically,  the 
best  move  from  h  randomly  generated  neighbors  is  selected. 

The  second  set  of  rules  diversify  the  search,  by  moving  the  local  search  into  a  new  region 
of  the  solution  space.  For  a  given  schedule  x,  the  new  region  is  reached  by  swapping  a  pair  of 
classes  of  the  same  course  between  two  pipes.  This  triple  (pipe  1,  pipe  2,  course)  is  the  best 
swap  among  m  randomly  generated  candidates.  To  prevent  the  search  from  cycling  back  to  the 
same  schedule  x,  once  a  swap  is  executed  it  will  be  tabu  for  the  next  n  iterations.  This  is  usually 
referred  to  as  intermediate  memory. 

The  third  set  of  rules  define  the  aspiration  criteria,  which  may  be  viewed  as  a  relaxation 
of  the  first  two  rules.  Since  the  potential  moves  using  the  above  rules  are  generated  randomly, 
it  is  possible  that  a  tabu  move  will  be  produced.  If  the  resulting  score  is  an  improvement  over 
the  best  schedule  found  so  far,  then  this  tabu  move  is  accepted. 


Given  a  schedule  x,  the  procedure  applies  the  first  set  of  rules  for  f  iterations.  This  is  fol¬ 
lowed  by  one  application  of  the  second  rule  (a  pipe  swap).  Since  an  infeasible  solution  may  be 
developed,  as  with  simulated  annealing,  a  very  large  penalty  is  assessed  for  violating  feasibility 
restrictions.  We  begin  the  tabu  search  with  the  first  feasible  solution  developed  by  the  greedy 
heuristic.  Additional  information  about  this  general  approach  may  be  found  in  Glover  [1989, 
1990a,  1990b]  and  Glover  and  Laguna  [1993]. 

2.2.4  The  Genetic  Algorithm 

The  genetic  algorithm  begins  with  a  schedule  generated  by  the  initial  application  of  the 
greedy  heuristic.  At  the  start  of  the  algorithm,  this  schedule  is  the  lone  member  of  a  pool  of 
schedules  that  will  be  built  up  and  modified  from  iteration  to  iteration.  Each  schedule  m  the 
pool  is  evaluated  to  determine  how  well  it  satisfies  the  problem  requirements  and  is  assigned  a 
score  based  upon  that  evaluation. 

The  algorithm  uses  two  basic  operations,  mutation  and  cross-over .  Given  a  schedule,  a 
mutation  replaces  a  random  number  of  randomly  selected  classes  with  randomly  selected 
classes  of  the  corresponding  course  and  shift.  A  cross-over  is  slightly  more  complicated  and 
involves  combining  two  schedules.  Given  two  schedules,  randomly  select  two  points  for  inter 
change.  The  two  points  are  selected  so  that  they  occur  at  pipeline  boundaries.  Place  pipes  in 
the  first  schedule  prior  to  the  first  interchange  point  and  pipes  in  the  second  schedule  between 
the  first  interchange  point  and  the  second  interchange  point  into  the  result  schedule.  Then  fill 
out  the  result  schedule  with  pipes  (if  any)  and  other  classes  from  the  first  schedule  after  the 
second  interchange  point. 

An  iteration  of  the  genetic  algorithm  proceeds  as  follows: 

1.  Randomly  select  a  given  proportion  of  the  pool  to  be  subjected  to  the  mutation  operation.  (Se¬ 
lected  schedules  remain  in  the  pool.)  Add  the  schedule  resulting  from  the  mutation  operation 
to  the  pool.  The  random  selection  of  schedules  to  mutate  is  biased  according  to  the  schedule 
score.  Better  schedules  are  more  likely  to  be  mutated. 


2.  Randomly  select  a  given  proportion  of  the  pool  to  be  subjected  to  the  cross-over  operation. 

(As  with  mutation,  selected  schedules  remain  in  the  pool.)  Add  the  schedule  resulting  from  the 
cross-over  operation  to  the  pool.  The  random  selection  of  schedules  to  cross-over  is  also 
biased  according  to  the  schedule  score.  Better  schedules  are  more  likely  to  be  selected  for  the 
cross-over  operation. 

3.  Reduce  the  pool  to  its  maximum  allowable  size  by  randomly  removing  schedules.  (This  is  not 
necessary  in  the  first  few  iterations.)  The  random  selection  of  schedules  to  remove  is  biased 
according  to  the  schedule  score.  Better  schedules  are  less  likely  to  be  removed,  and  the  best 
schedule  is  never  removed. 

In  any  particular  iteration,  schedules  may  be  generated  that  violate  one  or  more  of  the 
problem  restrictions.  In  addition  to  those  stated  in  the  section  on  simulated  annealing  (the 
offset  and  the  class  overlap  restrictions),  the  limit  on  the  number  of  classes  of  a  particular 
course  and  shift  maybe  violated.  Very  large  penalties  for  these  constraint  violations  are  added 
to  schedule  scores  of  the  infeasible  schedules.  These  infeasible  schedules  are  added  to  the 
pool  when  they  are  created;  however,  they  are  much  more  likely  to  be  removed  from  the  pool 
when  its  size  is  reduced  at  the  end  of  each  iteration. 

The  algorithm  repeats  the  process  for  a  specified  number  of  iterations  and  then  reports 
the  best  schedule  found.  Additional  information  about  this  general  approach  may  be  found  in 
Winston  [1992]. 

2.2.5  The  Evolutionary  Programming  Algorithm 

Evolutionary  programming  is  a  procedure  developed  by  D.  B.  Fogel  of  ORINCON  Cor¬ 
poration  which  also  uses  ideas  from  biology  to  produce  heuristic  algorithms  for  a  variety  of  dis¬ 
crete  optimization  problems.  For  example,  an  algorithm  for  the  traveling  salesman  problem 
can  be  found  in  Fogel  [1988];  a  technique  for  solving  linear  systems  may  be  found  in  Fogel  and 
Atmar  [1990];  and  a  technique  for  training  neural  networks  can  be  found  in  Fogel,  Fogel,  and 
Porto  [1990].  J.  P.  Lemoine  [1994]  adapted  evolutionary  programming  for  the  pipeline  sched¬ 
uling  problem  and  developed  a  software  implementation. 


Using  the  greedy  algorithm  to  obtain  starting  schedules,  Lemoine  develops  an  initial  set 
of  parent  schedules  which  are  then  perturbed  by  adding  random  gaps.  The  scores  for  each  par¬ 
ent  schedule  are  used  in  a  competition  based  upon  random  numbers,  that  determines  which 
parents  survive  to  become  the  basis  for  the  next  generation.  The  survivors  create  progeny  (a 
new  generation  of  schedules)  and  the  process  repeats.  The  progeny  are  produced  using  two 
methods.  One  routine  randomly  adds  and  removes  gaps  to  a  schedule.  The  other  randomly 
switches  positions  of  two  classes  in  a  randomly  selected  pipe. 

For  this  application,  a  generation  begins  with  20  schedules.  A  competition  reduces  these 
20  schedules  down  to  four  survivors.  Each  of  these  four  survivors  produces  five  progeny  result¬ 
ing  in  a  new  generation  of  20  schedules.  Lemoine  found  that  good  solutions  could  generally  be 
obtained  in  about  ten  generations. 

2.2.6  Empirical  Analysis 

Six  actual  problems  from  various  C— Schools  are  used  in  the  empirical  tests.  The  most 
difficult  problem  (PAS_2)  has  ten  courses  and  eight  pipelines  (see  Table  4).  All  problems  have 
a  one  year  planning  horizon  (50  weeks)  and  no  course  is  taught  by  more  than  three  shifts.  The 
course  length  varies  from  one  to  31  weeks. 

Since  Pipeline  Assistant  is  written  in  Turbo  Pascal  and  uses  routines  from  Microsoft  Win¬ 
dows,  the  first  test  was  conducted  on  a  486— based  PC.  A  comparison  of  Pipeline  Assistant  with 
a  C  language  implementation  of  the  greedy  algorithm  may  be  found  in  Table  4.  The  times  are 
wall  clock  times  and  are  rounded  to  the  nearest  second.  For  every  problem,  the  greedy  imple¬ 
mentation  produced  a  better  solution  in  less  time.  Since  Pipeline  Assistant  uses  a  random 
number  seed  based  upon  the  system  clock,  additional  runs  yield  different  solutions.  In  five 
additional  runs  with  FFG— 7,  Pipeline  Assistant  produced  schedules  having  the  following 
scores:  27, 34, 28, 26,  and  20. 


Table  4.  About  Here 


An  empirical  analysis  of  C  language  implementations  of  the  other  algorithms  may  be 
found  in  Table  5 .  Reported  times  are  user  times  (estimates  cpu  time)  on  a  DecStation  5000/240. 
For  each  of  the  interchange  heuristics,  the  initial  solution,  the  final  solution,  and  the  total  time 
is  reported.  Both  the  tabu  search  and  the  genetic  algorithm  apply  an  elementary  version  of  the 
greedy  method  to  obtain  initial  schedules.  The  evolutionary  programming  algorithm  applies  a 
more  advanced  version  of  the  greedy  algorithm  to  obtain  the  initial  schedule.  The  greedy  algo¬ 
rithm  produced  the  best  schedule  for  each  of  the  six  problems  and  was  always  faster  than  tabu 
search  and  the  genetic  algorithm.  The  evolutionary  programming  code  ran  each  problem  in 
approximately  17  seconds  regardless  of  the  problem  size.  On  PAS_2  the  evolutionary  pro¬ 
gramming  implementation  ran  faster  than  the  greedy  implementation,  but  failed  to  produce  as 
good  a  solution. 

The  evolutionary  programming  software  begins  by  generating  20  schedules  using  the 
greedy  method  with  local  improvement.  This  is  followed  by  15  iterations  of  the  evolutionary 
programming  algorithm.  The  score  of  the  best  of  the  greedy  schedules  is  reported  in  the  row 
entitled  initial  score.  Notice  that  the  computational  machinery  of  evolutionary  programming 
never  discovered  a  schedule  that  was  superior  to  that  produced  by  the  20  trials  of  the  greedy 
algorithm.  Hence,  the  best  solution  obtained  by  the  evolutionary  programming  software  was 
obtained  in  the  first  second  by  the  greedy  heuristic. 


Table  5.  About  Here 


2.3  An  Implicit  Enumeration  Algorithm  for  the  Pipeline  Scheduling  Problem 

In  this  Section  we  present  an  enumeration  algorithm  for  the  pipeline  scheduling  prob¬ 
lem.  A  lower  bound  is  developed  from  an  analysis  of  the  input  data  and  feasible  schedules 
(upper  bounds)  are  generated  one  pipe  at  a  time  using  a  specialized  enumeration  strategy. 
Hence,  a  partial  schedule  whose  score  is  greater  than  or  equal  to  the  incumbent  (best  schedule 
found  so  far)  can  be  fathomed  (discarded).  The  hope  is  that  a  good  incumbent  will  permit  the 
algorithm  to  fathom  most  of  the  possible  schedules. 

2.3.1  The  Lower  Bound 

Suppose  there  are  P  pipes  and  suppose  the  minimum  duration  of  each  of  these  pipes 
is  D  weeks.  That  is,  if  a  pipe  has  all  classes  scheduled  back— to— back,  then  the  pipe  lasts 
D  weeks.  Let  Tp  denote  the  target  start  date  for  pipe  p.  If  pipe  p  begins  on  or  after  week 
Tp  there  is  no  penalty.  Let  H  denote  the  planning  horizon  (either  50  or  100  weeks).  If  a 
pipe  lasts  through  week  H+t,  then  the  score  for  this  pipe  is  at  least  t.  Therefore,  L  = 
max[0,  H  -  (Ti  +  D  -  1)]  +  max[0,  H  -  (T2  +  D  -  1)]  +  ...  +  max[0,  H  -  (TP  +  D  - 
1)]  is  a  lower  bound  for  any  feasible  schedule. 

2.3.2  Complete  Enumeration 

In  order  to  develop  an  implicit  enumeration  algorithm,  one  needs  a  technique  to 
completely  enumerate  all  possible  solutions  since  in  the  worst  case,  the  algorithm  will  re¬ 
quire  a  complete  enumeration  of  the  solution  space.  Consider  a  problem  having  three 
courses  A,  B,  C  and  a  single  pipe.  There  are  3!  =6  possible  schedules  given  by  the  per¬ 
mutations  (ABC,  ACB,  BAC,  BCA,  CAB,  CBA).  That  is,  given  a  permutation,  convening 
dates  for  the  classes  are  determined  by  the  earliest  possible  time  for  each  class  in  the  per¬ 
mutation.  For  the  permutation  BAC,  class  B  is  convened  at  the  earliest  possible  week,  the 
target  week  for  the  pipe.  Class  A  begins  as  soon  as  B  is  completed  followed  by  C.  If  the 
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schedule  has  three  courses  and  two  pipes,  then  there  are  (3!)(3!)  =  36  possible  schedules  as 
illustrated  in  Figure  4.  In  general  for  a  schedule  having  C  courses  and  P  pipes  there  will  be 
(C!)p  possible  schedules  formed  in  this  way. 


Figure  4.  About  Here 


Note  that  the  tree  in  Figure  4  has  6+36=42  nodes.  The  first  six  nodes  represent  the 
first  pipe  and  corresponding  priorities  for  the  three  courses,  whereas  the  next  level  of  nodes 
represent  the  second  pipe  and  corresponding  priorities  for  its  courses.  In  general  a  tree 
will  have  (C!)  +  (C!)2  +  (C!)3  +  ...  +  (C!)p  nodes.  A  very  compact  algorithm  based  upon 
recursion  has  been  developed  to  generate  these  nodes.  A  node  is  defined  by  the  pair  (Z,p) 
where  p  denotes  the  pipe  number  and  Z  is  a  P  by  C  matrix  where  Zpj  denotes  the  course  in 
the  ith  position  of  pipe  p.  The  complete  enumeration  algorithm  is  given  as  follows: 
algorithm  complete  enumeration(P,  C) 

Inputs:  P  -  denotes  the  number  of  pipes 

C  -  denotes  the  number  of  courses 

Output:  A  set  of  (C!)p  matrices  Z  corresponding  to  all  of  the  permutations 

begin 

for  p  =  1  to  P  do  for  i  =  1  to  C  do  Zpj  -t—  i; 
perm(l,  1,  Z); 
end 
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procedure  perm(p,  i,  Z) 

Inputs:  p  —  denotes  the  current  pipe  • 

i  -  denotes  the  current  course  position 
Output:  A  set  of  Z  matrices 

begin  • 

if  i  <  C  then 

if  i  <  P  then  perm(p+l,  1,  Z); 

else  write  Z;  • 

else 

for  j  =  i  to  C  do 

switch(p,  i,  j,  Z);  • 

perm(p,  i,  Z); 
switch(p,  i,  j,  Z); 

end  ® 

end 
end 


procedure  switch(p,  i,  j,  Z) 

begin 

t  -t —  Zpj;  Zpi  -•  Zpjj  Zpj  ■<  t: 

end 


The  P  row  C  column  matrix  Z  will  correspond  to  one  permutation  of  the  P  pipes.  Hence,  the 
first  six  permutations  illustrated  in  Figure  4  correspond  to 


I"  1 2  3  1  fl23iri23l  r-1  fl23lfl23l 

L  123  J  L  132  J  L213  J  L23lJ  ^.312  J  L321  J 
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2.3.3  Construction  of  a  Pipe 

For  a  pipe  with  a  given  permutation  (i.  e.  Z),  the  convening  dates  for  the  classes  in 
the  pipe  are  determined  sequentially.  Given  a  partial  schedule  (i.  e.,  some  classes  have 
been  scheduled)  and  a  new  class  to  be  scheduled,  we  need  an  efficient  method  to  deter¬ 
mine  the  earliest  possible  convening  date  for  this  new  class.  Since  the  C— School  courses 
are  generally  team  taught,  they  frequently  have  multiple  offerings  of  the  same  course  being 
offered  simultaneously,  as  long  as  the  starting  dates  maintain  the  offset.  Suppose  course  A 
having  a  two  week  offset  currently  has  three  offerings  scheduled  as  illustrated  in  Figure  5. 
Then  weeks  8,  12, 15,  and  20  are  convening  dates  which  we  wish  to  consider.  That  is,  we 
consider  the  earliest  possible  starting  time  and  the  convening  times  plus  the  offset  for  the 
other  three  classes.  Let  y  denote  a  possible  trial  value  for  class  4.  If  y  is  in  the  open  inter¬ 
val  (10—2,10+2),  then  y  will  violate  the  offset  restriction  for  class  1.  Likewise,  if  y  is  in  the 
open  interval  (13—2, 13+2),  then  y  will  violate  the  offset  restriction  for  class  2.  If  Aj  is  the 
convening  time  for  class  i  with  an  offset  of  a,  then  y  can  not  be  in  the  open  interval  (Aj  — 
a,  Aj  +  a)  for  any  i.  These  ideas  result  in  the  following  assignment  algorithm: 
algorithm  assign(n,  A,  a,  (3,  y) 

Inputs:  n  —  denotes  the  number  of  currently  scheduled  classes  for  this  course 
Aj  —  denotes  the  convening  week  for  class  i 
a  —  denotes  the  offset 

(3  —  denotes  an  earliest  possible  convening  week 
Output:  y  —  convening  week  for  class  n+1 
Assumption:  Ai  <.  A2  <L ...  <L  An 

begin 

y-^P; 

for  i  =  1  to  n  do  if  A,  -  a  <  y  <  Aj  +  a  then  y  -t—  max(y,  Aj  +  a); 

end 
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Given  a  permutation  of  pipes  and  classes  defined  by  Z,  assign  can  be  used  to  construct  a 
schedule.  That  is,  assign  determines  the  convening  time  for  each  class  in  Z. 

2.3.4  Implicit  Enumeration  Algorithm 

Given  a  permutation  matrix  Z  and  a  p  <  P,  a  partial  schedule  can  be  developed  by 
applying  assign  sequentially  to  the  first  p  pipes.  Let  score(Z,p)  denote  the  score  of  the  par¬ 
tial  schedule  defined  by  Z  and  p.  Since  score(Z,p)  <.  score(Z,p+l)  then  score(Z,p)  is  a 
lower  bound  on  score(Z,P).  Combining  this  notion  with  the  complete  enumeration  algo¬ 
rithm  results  in  the  following  implicit  enumeration  algorithm: 
algorithm  implicit  enumeration(P,  C) 

Inputs:  P  —  denotes  the  number  of  pipes 

C  —  denotes  the  number  of  courses 
Output:  v  —  denotes  the  score  of  the  best  solution  found 
begin 

for  p  =  1  to  P  do  for  i  =  1  to  C  do  Zpj  +—  i; 

let  v  denote  the  score  of  the  solution  found  using  the  greedy  heuristic; 

let  L  denote  a  lower  bound  based  on  the  analysis  given  in  Section  2.3.1; 

if  v  >  L  then  enumerate(l,  1,  Z,  v); 

write  “the  best  schedule  found  has  a  score  of  v” 


procedure  enumerate  (p,  i,  Z,  v) 

Inputs:  p  —  denotes  the  current  pipe 

i  —  denotes  the  current  course  position 
Z  -  denotes  a  permutation 
v  —  denotes  the  current  best  solution 

begin 

if  i  <  C  then 
if  p  <  P  then 

if  score(Z,p)  <  v  then  enumerate(p+l,  1,  Z,  v); 
else  if  score(Z,P)  <  v  then  v  -t—  score(Z,P); 

end 

else 

for  j  =  i  to  C  do 
switch(p,  i,  j,  Z); 
enumerate(p,  i,  Z,  v); 
switch(p,  i,  j,  Z); 

end 

end 

end 

Of  course,  one  must  also  save  the  schedule  whenever  v  is  updated.  The  hope  is  that  the  initial 
solution  (upper  bound)  will  be  small  and  that  the  lower  bounds  will  be  large. 

2.3.5  Empirical  Analysis 

The  six  test  problems  were  run  using  a  C  language  implementation  on  a  275  MHz  Alpha 
based  DecStation  and  the  results  are  summarized  in  Table  6.  For  the  first  four  problems,  the 
greedy  heuristic  found  a  solution  whose  value  was  equal  to  the  lower  bound;  hence,  no  enumer¬ 
ation  was  needed.  A  provable  optimum  for  the  FFG— 7  problem  was  obtained  after  approxi- 


mately  one -half  hour  of  computer  time.  On  PAS_2  the  system  was  terminated  after  it  ex¬ 
amined  three  million  nodes  which  took  approximately  one  hour  of  cpu  time.  At  termination, 
no  improvement  over  the  greedy  heuristic  had  been  made  and  the  gap  between  the  upper  and 
lower  bounds  was  approximately  25%. 


Table  6.  About  Here 


2.3.6  Recommendation 

Based  on  our  empirical  analysis,  we  recommend  an  algorithm  that  uses  three  proce¬ 
dures:  (i)  the  greedy  heuristic,  (ii)  the  lower  bound,  and  (iii)  the  implicit  enumeration  algo¬ 
rithm.  The  greedy  heuristic  has  consistently  produced  fairly  good  schedules  which  could  be 
used  by  the  client.  The  lower  bound  should  also  be  produced  and  reported  to  the  user.  If  the 
gap  between  the  upper  and  lower  bound  is  large,  then  the  user  should  be  given  the  option  to  run 
the  implicit  enumeration  algorithm  for  some  given  time,  to  determine  if  the  gap  can  be  closed. 
A  finite  termination  criteria  should  be  incorporated  to  prevent  excessive  run  time  in  the  im¬ 
plicit  enumeration  algorithm. 


III.  NAVY  A- SCHOOLS 


The  basic  courses  for  new  recruits  are  offered  at  campuses  known  as  A-Schools.  The 
largest  are  located  at  Norfolk,  San  Diego,  and  Great  Lakes,  and  during  FY96  over  80,000  sail¬ 
ors  took  at  least  one  A-School  course.  The  managers  of  each  A— School  is  responsible  for 
planning  the  class  schedule  for  their  School.  In  this  Section  we  describe  the  A-School  class 
scheduling  problem,  we  present  a  greedy  heuristic  for  this  problem,  and  we  provide  an  empiri¬ 
cal  evaluation  of  our  algorithm. 

3.1  Description  of  the  Problem 

A  Navy  A— School  will  offer  a  variety  of  training  courses  during  a  one  year  planning  peri¬ 
od.  Some  courses  last  only  a  few  days  while  others  may  involve  over  one  hundred  days  of 
instruction.  Some  days  a  class  may  meet  in  a  large  lecture  hall  and  other  days  the  class  may  be 
divided  into  smaller  groups.  These  smaller  groups  will  meet  in  different  rooms  with  individual 
instructors.  Hence,  on  day  one,  the  School  may  need  a  large  lecture  hall  and  one  instructor  but 
on  day  two,  the  School  may  need  three  small  lecture  rooms  and  three  instructors.  The  instruc¬ 
tors  are  interchangeable  so  that  an  instructor  can  teach  either  the  whole  group  in  a  large  lec¬ 
ture  hall  or  the  smaller  groups.  In  addition,  some  classes  may  require  special  equipment  such 
as  calculators  or  personal  computers.  The  managers  of  the  A-  Schools  are  responsible  for  de¬ 
veloping  a  class  schedule  that  requires  the  minimum  number  of  instructors  subject  to  resource 
constraints  on  lecture  halls,  seminar  rooms,  laboratory  space,  and  equipment. 

The  input  data  used  to  describe  the  A-School  class  scheduling  problem  is  given  below: 
Let  i  -  denote  a  course  (i=  l,...,n) 

Let  Dj  —  denote  the  duration  of  course  i  in  days.  (This  can  vary  from  a  few  days  to  well  over 
one  hundred  days.) 

Let  Oj  -  denote  the  number  of  times  course  i  will  be  offered  during  the  year.  All  offerings 


could  convene  on  the  same  day  or  they  could  be  dispersed  throughout  the  year .  Dis 
persal  usually  results  in  the  requirement  of  fewer  instructors  on  the  busiest  day. 

Let  Ijd  —  denote  the  number  of  instructors  needed  by  course  i  on  day  d.  This  can  vary  from 

a  minimum  of  one  to  a  maximum  of  four. 

Let  Cjid  -  denote  the  consumption  of  resource  j  by  course  i  on  day  d.  These  are  usually  0  or  1 

depending  upon  whether  a  particular  type  room  is  needed. 

Let  Rjt  —  denote  the  number  of  resources  of  type  j  available  on  calendar  day  t.  Usually  the 

consumption  is  a  small  integer  corresponding  to  the  number  of  rooms  of  a  given 
type  available  at  the  School. 

The  decision  variables  are  convening  dates  for  the  Oi  +  O2  +  ...  +  On  classes  to  be  of¬ 
fered.  The  objective  is  to  select  convening  dates  that  require  the  least  number  of  instructors  on 
the  busiest  day  of  the  year.  The  problem  may  be  further  complicated  by  the  imposition  of  con¬ 
vening  day  restrictions.  For  example,  some  A— Schools  require  that  a  certain  course  always 
begin  on  a  Monday.  Some  may  restrict  the  number  of  offerings  of  a  given  course  in  a  given 
month. 

3.2  The  Greedy  Algorithm 

Consider  an  A- School  which  offers  Oi  +  02  +  ...  +  On  =  m  classes  during  the  year. 
Suppose  this  School  is  open  Monday  through  Friday  except  for  thirteen  holidays.  For  a  non- 
leap  year,  this  results  in  356  -  104  -  13  =  248  possible  convening  days  for  each  of  m  classes. 
Hence,  there  are  248m  possible  schedules. 

Based  upon  the  successful  implementation  of  the  greedy  algorithm  for  C-Schools,  we 
developed  a  similar  approach  for  the  A— School  scheduling  problem.  Consider  a  School  that 
offers  three  courses,  A,  B,  and  C.  We  first  develop  the  6  permutations  for  these  courses  (ABC, 
ACB,  BAC,  BCA,  CAB,  CBA).  Suppose  that  A  is  offered  twice,  B  once,  and  C  three  times. 
Using  a  simple  rule  to  spread  the  offerings  of  each  course  throughout  the  year,  we  obtain  the 
following  possible  sequences: 
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Sequence  Number 
1 


Class  Order 


Al— B1-C1-C2-A2-C3 

2  A1-A2-C1-C2-B1-C3 

3  B1-A1-C1-A2-C2-C3 

4  B1-C1-A1-C2-A2-C3 

5  Cl— A1-B1-C2-A2-C3 

6  Cl— B1-A1-C2-A2-C3 

For  sequence  1,  the  greedy  algorithm  will  schedule  the  first  offering  of  A  followed  by  the  first 
offering  of  B.  This  is  followed  by  two  offerings  of  course  C,  the  second  offering  of  A,  and  the 
final  offering  of  C. 

The  greedy  algorithm  initializes  the  day  index  t  to  1  and  assigns  the  first  class  in  the  se¬ 
quence  a  convening  day  of  t.  For  sequence  1,  Al  convenes  at  day  1 .  Next  it  attempts  to  assign 
the  second  class  in  the  sequence  a  convening  day  of  t.  If  a  resource  violation  occurs,  then  t  is 
incremented  by  one  and  it  tries  again.  Once  the  second  class  is  scheduled,  then  it  proceeds  to 
the  third.  This  process  continues  until  either  all  classes  are  scheduled  within  the  year  (a  feasi¬ 
ble  schedule)  or  it  fails  to  find  a  feasible  schedule.  If  a  feasible  schedule  is  obtained,  then  an 
upper  bound  on  the  number  of  instructors  needed  has  been  established.  Suppose  this  process 
results  in  a  feasible  schedule  using  50  instructors  on  the  busiest  day.  This  first  schedule  found 
does  not  attempt  to  minimize  the  number  of  instructors  and  the  algorithm  is  concerned  only 
with  feasibility.  In  the  next  stage  the  number  of  instructors  is  systematically  reduced  using  a 
binary  search  on  the  interval  [1,50]  and  the  algorithm  proceeds  to  find  a  feasible  schedule  with 
this  sequence  and  the  minimum  number  of  instructors.  This  strategy  is  repeated  with  the  rest 
of  the  sequences  and  the  best  of  the  best  is  saved. 


3.2  Empirical  Analysis 


Seven  problems  were  used  in  our  empirical  tests.  The  number  of  courses  varies  from  one 
to  six  and  the  number  of  resource  restrictions  varies  from  one  to  nine.  The  course  length  varies 
from  ten  days  to  over  one  hundred  days  and  the  number  of  offerings  of  a  single  course  during 
the  year  varies  from  one  to  sixty.  The  characteristics  of  these  test  problems  may  be  found  in 
Table  7. 


Table  7.  About  Here 


SABRE  Decision  Technologies  has  developed  a  PC  based  system  to  obtain  schedules  for 
this  problem.  The  system  is  called  NCSS  (Navy  Class  Scheduling  System)  and  uses  simulated 
annealing  in  an  attempt  to  obtain  a  good  schedule.  Table  8  summarizes  our  results  using  both 
simulated  annealing  and  the  greedy  heuristic  on  the  seven  test  problems.  On  six  of  the  prob¬ 
lems  ,  greedy  provided  a  better  answer  and  on  the  seventh  there  was  a  tie.  Since  NCSS  is 
written  in  Visual  Basic,  it  was  run  on  a  50  MHz  486  based  PC  and  each  run  took  approximately 
15  minutes.  The  greedy  code  is  written  in  C  and  those  runs  were  made  on  a  275  MHz  Alpha 
based  DecStation  and  never  required  more  than  a  few  seconds..  These  results  are  consistent 
with  those  obtained  previously  with  the  pipeline  scheduling  problem. 


Table  8.  About  Here 


IV.  SUMMARY  AND  CONCLUSIONS 


Five  computer  codes  designed  to  obtain  good  solutions  to  the  Navy’s  C— School  schedul¬ 
ing  problem  involving  pipelines  have  been  tested.  The  five  codes  are  based  upon  the  following 
algorithms: 

(i)  simulated  annealing, 

(ii)  a  greedy  heuristic, 

(iii)  tabu  search, 

(iv)  a  genetic  algorithm,  and 

(v)  evolutionary  programming. 

In  empirical  tests  with  these  five  codes  on  six  actual  C— School  scheduling  problems,  the 
greedy  algorithm  produced  the  best  schedules.  Not  only  did  it  produce  the  best  schedules,  its 
computational  time  was  least  for  all  runs  except  for  problem  PAS_2.  For  this  problem,  the  evo¬ 
lutionary  programming  code  terminated  after  18  seconds  while  the  greedy  algorithm  ran  for 
about  one  minute. 

The  four  interchange  heuristics  (simulated  annealing,  tabu  search,  the  genetic  algo¬ 
rithm,  and  evolutionary  programming)  all  involve  several  tuning  parameters  (initial  tempera¬ 
ture,  iterations  before  a  temperature  change,  number  of  iterations  before  the  local  search  is 
abandoned  and  a  new  region  is  explored,  the  mutation  and  cross-over  rules  used,  the  number 
of  schedules  retained  in  a  generation,  the  number  of  generations  computed,  et  cetera).  These 
tuning  parameters  were  set  by  the  authors  of  each  code  and  were  not  modified  in  our  experi¬ 
ments.  However,  it  is  well  known  that  modification  of  the  tuning  parameters  in  any  of  these 
codes  may  produce  different  results.  The  only  tuning  parameter  in  the  greedy  algorithm  is  the 
maximum  number  of  schedules  to  be  developed  before  termination.  This  parameter  which 
was  set  to  10,000  in  these  runs  could  probably  be  set  to  a  much  smaller  value  with  no  deleteri- 


ous  effect  on  the  final  solution.  In  the  empirical  analysis,  the  best  solution  for  the  problem 
PAS_2  was  discovered  at  iteration  1746.  It  is  also  possible  that  all  of  the  interchange  heuristics 
could  benefit  by  first  running  the  greedy  heuristic  to  completion  followed  by  the  interchange 
search. 

We  also  developed  an  implicit  enumeration  algorithm  for  the  pipeline  scheduling  prob¬ 
lem.  At  the  expense  of  substantial  computational  time,  this  algorithm  discovered  the  optimal 
solution  for  the  FFG-7  problem.  None  of  the  heuristic  methods  were  able  to  find  this  optimal 
solution. 

The  ideas  used  in  the  successful  greedy  heuristic  for  the  pipeline  scheduling  problem 
were  extended  for  the  A- School  scheduling  problem.  In  empirical  tests,  the  greedy  heuristic 
performed  better  than  an  implementation  of  simulated  annealing  on  a  set  of  test  problems. 
Based  on  these  studies,  we  conclude  that  greedy  type  algorithms  are  ideally  suited  for  class 
scheduling  problem  for  Navy  training  schools. 
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Weeks 


Figure  1.  Schedule  for  a  Ten  Week  Course  Having  a  TWo  Week  Offset 
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Figure  2.  Example  of  TWo  Pipes 
(Pipe  1  =  [ABC]  has  no  gaps) 

(Pipe  2  =  [CAJB]  has  four  weeks  of  gaps) 
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Weeks 
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Figure  5.  Potential  Convening  Dates  for  a  Fourth  Offering  of  Course  A 
(Earliest  Possible  Convening  Date  is  Week  8  and  the  Offset  is  2) 


Table  1.  The  FFG-7  Pipeline  Scheduling  Problem  (Actual  Input  Data  File) 
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Table  3.  Pipeline  Scores  for  the  Schedule  Illustrated  in  Figure  3 
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Table  4.  Computational  Comparison  of  Pipeline  Assistant  with  the  Greedy  Algorithm 

(All  times  are  in  seconds  on  a  50  MHz  486 -based  PC.) 


B-38 


*No  feasible  solution  found. 


Table  5.  Empirical  Analysis  of  the  Algorithms 
(All  times  are  in  seconds  on  a  DecStation  5000/240.) 
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Evolutionary  Programming 

initial  score  0  0  3  28  21  293 

final  score  0  0  3  28  21  293 

total  time  18  18  16  16  17  18 


Table  6.  Empirical  Analysis  of  the  Implicit  Enumeration  Algorithm 
(All  times  are  in  seconds  on  an  Alpha  based  DecStation.) 
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^terminated  after  3,000,000  nodes 


Table  7.  Characteristics  of  Test  Problems 


Table  8.  Comparison  of  the  Simulated  Annealing  Algorithm  and 
the  Greedy  Heuristic  on  the  A- School  Scheduling  Problem 
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Generalized  Networks 
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ABSTRACT 

This  Chapter  describes  a  specialization  of  the  dual  simplex  algorithm  for  the  gener¬ 
alized  network  problem.  We  use  a  dual  two  phase  method  along  with  efficient  dual 
partial  pricing  schemes  and  specialized  routines  for  the  dual  ratio  test.  In  comparison 
with  CPLEX  3.0  on  a  set  of  ten  benchmark  problems,  we  found  that  our  special¬ 
ized  dual  code  performed  20  times  faster  than  the  best  CPLEX  optimizer.  Problems 
having  10  to  20  thousand  arcs  are  routinely  solved  in  under  10  seconds  on  a  60MHz 
DECStation  5000/260  with  our  rapid  dual  generalized  network  optimizer. 


1  INTRODUCTION 

The  mathematical  problem  of  finding  a  minimal  cost  flow  through  a  capacitated 
network  is  a  fundamental  problem  in  both  operations  research  and  computer 
science.  A  generalization  of  this  model  which  allows  for  either  gains  or  losses 
as  flows  pass  along  an  arc  is  known  as  the  generalized  network  problem.  Gen¬ 
eralized  networks  have  been  used  to  model  a  wide  array  of  applications  in 
many  engineering  and  economic  areas.  Applications  that  involve  either  gains 
or  losses  include  electric  power  carried  on  transmission  lines,  cash  management 
that  involves  the  time  value  of  money,  and  manufacturing  processes  of  vary¬ 
ing  efficiencies.  Other  applications  allow  for  flow  conversion  from  one  unit  to 
another.  Further  discussion  of  a  wide  variety  of  applications  can  be  found  in 
Ahuja,  Magnanti,  and  Orlin  [1]  ,  Glover,  Klingman,  and  Phillips  [11],  Glover 
et  al  [10],  and  Mulvey  and  Zenios  [19]. 
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1.1  Problem  Description 

Let  efc  -  (i  j)  denote  an  arc  from  node  i  to  node  j  having  multiplier  ak  associ¬ 
ated  with  node  i,  and  6,  associated  with  node;  If  i  =  j,  then  arc  efc  has  one  end 
which  is  UM  on  nod,  i  with  yi^-t,phe,  « ^  V  -  (U, .. 

1^<c»'h-.'Kp^donadh,et,d6,.ph 

<  V}A  >  For  ejp  =  (*,  j)  let 


k 

9i 


ak ,  for  t  —  i\ 
bk,  for  f  =  ;; 
0,  otherwise. 


Tptr_  rfli,02|  \an],  be  an  m  x  n  matrix  called  the  node-arc  coefficient  ma- 

bounds  ,arc  capacities,  and  arc  >  P  -  _  r  _  spn  .  qx  —  r  l  <  x  <  u} 

generalized  network  problem  can  be 
stated  mathematically  as  min*  {cx  :  x  €  $}• 


1.2  Survey  of  Literature 


Generalized  networks  ha»e  the  special  structure  that  the  basis  can  be ^pre- 

K,nted  r  “Lt:  SiSSU  «* 

(NETOMhat  exploits  the  ^  GtovTc' ri  a  [10]“nh“S 

SO  for  the^starting  strategy,  the  pivot  s, 

fertfon  criteria,  and  degeneracy  handling.  In  1984,  Brown  and  McBr.de  [6 
printed  a  detailed  description  of  a  complete  implementation  of  “  '®c“" 
nrimal  simplex  system  specialized  for  the  generalized  network  problem.  Th 
s?Zm  „  2s  a  preorder  traversal  method  in  addition  to  the  predecessor,  depth 

nWkated  the  performance  efficiency  of  the  primal  simplex  generalized  net¬ 
work  code  LPNETG,  when  using  different  internal  programming  tactics  sue 
no  altram ative  nivot  strategies,  column  normalization,  and  the  big-M  starting 
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method.  In  1988,  Nulty  and  Trick  [22]  developed  the  first  primal  simplex  code 
written  in  the  C  language.  They  use  the  predecessor,  thread,  reverse  thread, 
and  the  level  node  labels  presented  in  Kennington  and  Helgason  [15]  to  repre¬ 
sent  the  basis.  In  1988,  Bertsekas  and  Tseng  [5]  proposed  a  new  class  of  al¬ 
gorithms,  that  adopt  the  nonlinear  programming  relaxation  methods,  to  solve 
the  linear  cost  generalized  network  problem.  The  algorithm,  is  based  on  the 
iterative  improvement  of  the  dual  cost  while  maintaining  a  flow  that  satis  es 
complementary  slackness.  In  1992,  Clark  ei  al.  [7]  developed  and  empirically 
tested  two  primal  simplex  parallel  algorithms,  GENFLO,  and  GRNET2,  for 
solving  generalized  networks.  Both  codes  use  a  gradient  penalty  method  to 
find  a  starting  feasible  solution. 


1.3  Objective  of  Investigation 

Although  several  papers  have  been  presented  that  discuss  the  development 
of  specialized  algorithms  and  software  implementations  for  solving  generalized 
networks,  none  of  these  are  based  on  the  dual  simplex  method.  The  objective 
of  this  investigation  is  to  develop  and  perform  an  empirical  analysis  of  a  spe¬ 
cialization  of  a  dual  simplex  algorithm  for  the  generalized  network  problem. 
There  are  two  key  issues  in  the  development  of  a  dual  procedure:  (l)  the  cre¬ 
ation  of  an  effective  dual  pricing  scheme,  and  (ii)  the  creation  of  an  efficient 
scheme  for  computing  an  updated  row  of  the  coefficient  matrix  which  exploits 
the  underlying  network  structure. 


2  THE  DUAL  SIMPLEX  ALGORITHM 

Duality  theory  in  linear  programming  is  well  known.  Each  variable  in  the  pri¬ 
mal  corresponds  to  a  constraint  in  the  dual,  and  each  constraint  in  the  prima  is 
associated  with  a  dual  variable.  The  dual  simplex  method  has  received  limited 
attention,  and  its  role  has  been  limited  to  sensitivity  analysis,  parametric  pro¬ 
gramming,  and  the  solution  of  integer  programming  problems  ( see  Nemhauser 
and  Wolsey  [21]  and  Parker  and  Rardin  [23]). 

There  are  some  instances  where  the  dual  method  has  an  advantage  over  the 
primal  method,  see  Ali,  Padman,  and  Thiagarajan  [2].  The  reported  results 
of  their  dual  simplex  code,  specialized  for  pure  network  problems,  illustrates 
improved  performance  on  a  subset  of  the  benchmark  of  standard  NETGEN 
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t  •  _  j  qtutz  [161),  and  an  additional  set  of  larger 

problems  (Klingman,  Napier,  and  Stuta  [  J h 

transportation  problems. 

Under  the  non-degeneracy  assumption^  ^  otf  with  an  im- 

with  a  basic  primal  Seated  until  primal  optimality  is 

proved  objective  functron.  These  pivot ^  dual  method  starts  with  a 
achieved,  i.e.  obtaining  a  dua  Qne  wit^  an  improved  objective  func 

dual  feasible  solution  and  Piv°ts  t0  .  ,  infeasibilities.  This  is  repeated  until 

tion,  while  attempting  to  reduce ^  p  Jn  this  section  we  present  the 

dual  optimality  (  primal  feasibil  y  )  f  selecting  the  leaving  variable, 

Let  *,  denote  the  dual 

with  the  upper  bound  constra  3  denote  the  set  of  dual  feasible 

3 i",M  €  *")  ••  G?r  +  A  '  J  ",  nrihTem  forlhe  generalized  network  problem  is 
solutions.  Therefore,  the  dual  problem  8 .  £  y>  let  q  =  A(q)  denote  the 

max^.x.M  {r»  +  lX~  ^  •{*>  and  let  e*  denote  the  corresponding 

index  of  the  Recall  that  G  is  an  m  x  n  matrix  with 

column  of  an  ro  x  m identi tty  matn: x  ^  ^  ^  denotes  a  basis.  Let 

full  row  rank,  g3  €  5ft  deno  variables  be  partitioned  into  basic  and 

x  be  an  extreme  point  and  let  the ^  one  ofit’s  bounds  i.e. 

ZXo  LetM  =  {*  :  *  =  ^ 

of  the  non-basic  variables^ ^^n-Wc’ variables  at  upper  bound,  and  let 

xN  e  denote  thbleVseClnd  tetT~=  ^W^enote'  the 

XB  e  denote  the  vector  of  b3310  varia  ’  €  sftm  denotes  the  vector  of 

set  of  indices  of  the  basic  variables.  Recall  «**<=*  ^  =  r  _  and  that  A(i) 
duals  associated  with  the  flow  c°“sf^a  ^  ^  associated  with  constraint  i. 

~  -  The  d"ex 

algorithm  can  be  stated  as  follows: 

Procedure  Dual 

Input :  c}  G,  UyB^Mu  «• 

Output :  xB ,  x^  j  v  . 

bee^-cBB-'\ n  <-«-*/,  ll=Xi 
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XiB  . 


■  EfceJV.  CkUk  +  Ylketf,  Cklk  +  CkXk 


v  _ _ _ 

I  <-  {i  £  V  :  (xf  <  If)  or  (xf  >  uf )}; 

while(  X  <j>  )do 

select  gel;  q  *-  A(q);  /  *  leaving  variable  *  / 

z  <_  e$ jj-i;  /  *  row  g  of  basis  inverse  *  / 

to  +-  sign(u?  -  xq)zG]  /  *  updated  row  *  / 

£  <_  {j  e  .A/)  ;  Wj  <  0}  U  {j  €  jV'u  :  wj  >  0}; 
if  (£  =  4>  )then  exit(‘unbounded  dual  problem’); 

maxjgf 


{  }’ 


y  * 
A 


xf 


argmaxjg£ 

V;  V  -  {i  :  Vi  #  0}; 

if  xq  <  lq  1 


f  iSszdsl  j 

j  i 

b  y 4 

—  xf  -  At/i,  i 


if  xq  >  uq 
-  At/i,  i  €  P;  x5 


■  x,  +  A; 


if(  S  <  0  )then 

cr*  < —  crfc  —  6tOfc,  k  £//]  Uq  *~ 
■  fi  +  6z;  v*  *—  v*  +  <r, A; 


6 

-<5 


if  Xg  = 
if  x„ 


/  *  ratio  test  *  / 
/  *  updated  column  *  / 

/  *  flow  update  *  / 

=*. 

=  Ug  J  ’ 


7T  < 

endif 

B  <—  5U{s}\{g};  JV)4— M\{«};  Mi 
if(  =  /g  )then  {g}; 

iff  )then  Mi  <—  Mi  U  {q}\ 

X  *-  (i  e  V  :  (xf  <  /f )  or  (xf  >  uf )}; 
endwhile 
end. 


/  *  duals  and  objective  update  *  / 
jVu\{s};  /  *  basis  update  *  / 


2.1  A  Specialized  Implementation  for 
Generalized  Networks 

The  main  operations  of  a  dual  simplex  pivot  are  selecting  the  leaving  variable, 
calculating  the  updated  row,  determining  the  entering  variable,  and  updating 
the  flows,  the  reduced  cost,  and  the  inverse  of  the  basis.  For  generalized  net¬ 
works,  most  of  these  operations  can  be  carried  out  on  the  tree  representing  the 
basis,  and  the  inverse  of  the  basis  is  never  calculated. 

Calculating  a  row  of  the  basis  inverse,  z  =  can  be  considered  as  a 

special  case  of  the  dual  calculations  where  the  vector  c  is  replaced  by  e*  with 
a  single  nonzero  entry  in  position  q.  For  more  details  on  the  specialized  dual 
calculations  see  Mohamed  [18].  Let  Z  —  {i :  z*  ^  0},  denote  the  set  of  indices 
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for  the  nonzero  entries  in  z.  Let  Bq  be  the  basis  partition  corresponding  to  the 
augmented  forest  component  that  contains  node  q ,  and  let  zq  and  e9,  be  the 
matching  partitions  of  z  and  e9.  Suppose  that  Bq  corresponds  to  a  rooted-tree, 
Tr.  The  system  zqBq  =  e9  can  be  solved  using  back  substitution,  and  Z  will 
only  contain  the  nodes  in  the  subtree  component  rooted  at  node  q.  Suppose 
that  Bq  corresponds  to  a  one-tree,  T°  and  consider  only  the  cycle  of  T° .  For 
notation  purpose,  assume  that  the  cycle  nodes  are  numbered  1,  First, 

the  zq  values  for  the  cycle  nodes  are  obtained  using 

- - - r,  and 

ai(l  -  Wi W2  •  *  ‘U>m) 

-CLi  •  ! 

— — Ziy  2  —  —  1. 

bi 

Then,  the  zq  values  for  the  remaining  nodes  in  the  one-tree  are  calculated  using 
back  substitution.  In  this  case  Z  will  contain  all  nodes  in  the  one-tree  T° . 

The  calculation  of  the  the  updated  row,  w  =  zG,  is  a  dominant  operation  in  the 
dual  pivot.  Preliminary  profile  investigation  revealed  that,  when  using  standard 
matrix  operations,  approximately  50%  of  the  time  is  spent  in  calculating  w  with 
a  complexity  of  0(mn).  This  time  can  be  drastically  reduced  by  exploiting  the 
sparsity  of  z  and  the  special  structure  of  the  network  problems.  Let  G*  denote 
the  iih  row  of  G.  To  exploit  the  sparsity  of  z,  the  updated  row  calculation  is 
carried  out  as  follows: 

w  =  J2ziGi. 

iez 

This  reduces  the  complexity  of  calculating  w  to  0(\Z\n),  where  \Z\  <  m  in 
most  cases.  Let  T  =  {k  :  ek  =  (t,  v)  G  Ay  v  G  V},  the  forward  star  of  node 
2,  and  let  T*  =  {k  :  e*  =  (u,  f)  €  -4,  v  G  V},  the  backward  star  of  node 
i.  Let  W  ~  {j  :  Wj  ^  0}.  For  network  problems,  G*  can  be  represented 
using  the  forward  and  the  backward  stars  of  node  i.  This  special  structure  of 
the  network  problem  has  a  significant  influence  on  the  sparsity  of  w}  that  is 
W  =  U izziF'  UP).  This  reduces  the  complexity  of  calculating  the  updated 
row  to  0(\Z\h)y  where  h  is  the  average  number  of  arcs  incident  to  a  node.  In 
practice  n<n. 

The  entering  variable,  xs ,  and  the  change  in  the  reduced  cost,  (5,  are  determined 
in  the  ratio  test,  which  is  conducted  on  the  nonbasic  variables  in  £  C  W.  The 
updated  column,  y}  is  calculated  using  the  specialized  procedure  described  in 
Helgason  and  Kennington  [15].  The  flow  update  is  carried  out  only  on  the 
basic  arcs  associated  with  the  nodes  in  V.  For  additional  information  see 
Mohamed  [18].  This  is  followed  by  updating  the  reduced  cost  of  the  nonbasic 
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variables  in  W,  updating  the  node  labels,  updating  B  and  N  to  obtain  the  new 
basis,  and  updating  the  duals. 


Let  k  denote  the  kth  iteration  of  the  dual  simplex  algorithm.  Let  be  the 
vector  of  duals  at  iteration  fc,  let  be  the  basis  matrix  at  iteration  and 
let  D W  =  [jBW]"1.  Let  q  E  B  be  the  leaving  variable,  and  let  q  denote  it’s 
position  in  B.  Therefore,  the  dual  vector  at  iteration  fc  +  1  is  defined  as  follows: 


T(fc+1)  =  cB'fc+1>  £,(*+1)  t 

(6.1) 

where 

£(*+1)  _  |j  _  J_(y  _  e?VT| 

(6.2) 

and  „  v 

(6.3) 

Substituting  (6.2)  and  (6.3)  into  (6.1)  we  obtain 


*(*+!)  =  j>fc)  -  cbW  (^2_ii)  e?‘e?‘T  )  | £><*>  -  i(2 /  -  e«)e^  £><*>) 

_  cB(l)  D(fc)  _  cB(fc)  f  c?  -  c A  eiefDW  _  cB(k)  i>y~eileirD^ 

\  Cq  J  Vq 

+C*W  (c^z£f.)  e^L {y-ei)eiTDW 
\  Cq  J  y? 

+cB“)  (^T1^1)  e«VT(j/-e«))e?'T£»(i) 

=  .‘  +  ^{cBW9t(^)e‘  +  c»'-’ («*-„) 

+cB(k)  (£t^fl)  -  1))  e*T£>(*> 

=  ^  +  1  j>W  (e*  —  y)  —  cs<fc) e«  (^i)  }  ^ 

=  ,‘ +  e«TD<‘) 
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=  *-«  +  (t-)  e?T£>(fc) 

Therefore,  the  dual  vector  can  be  updated  using  +  (j^j  zi  an<i 

the  sparsity  of  z  can  be  exploited  in  this  update. 


2.2  Dual  Pricing  Strategies 

Since  all  nonbasic  variables  assume  one  of  their  bounds,  primal  infeasibility  oc¬ 
curs  only  when  the  basic  variables  violate  the  bounds.  The  dual  simplex  pricing 
strategy  to  select  the  leaving  variable  xq  has  received  very  little  attention  in 
the  literature.  Although  the  selection  of  the  leaving  variable  appears  to  be  a 
simple  operation,  the  exact  strategy  used  is  very  important  in  the  development 
of  an  efficient  generalized  network  code. 

The  standard  approach  for  choosing  the  leaving  variable  is  to  scan  the  basic 
variables  and  determine  the  variable  having  maximum  bound  violation.  The 
index  of  the  leaving  variable  selected  using  the  maximum-infeasibiliiy  strategy 
is 

J  h-Xk  if  Xk  <  k  1 
5  =  argmaxfe6B|  if  }  ■ 

Using  this  approach,  every  basic  variable  is  scanned  at  each  iteration,  which 
requires  a  node  length  search.  This  section  presents  techniques  which  improve 
this  search  by  performing  a  partial  pricing  where  only  a  subset  of  the  basic 
variables  are  scanned. 

Let  p  <  |V|  denote  the  page  length.  The  fixed-page  pricing  technique  scans  a 
page  of  basic  variables  at  a  time.  If  one  prices  favorable,  then  it  is  selected; 
otherwise,  another  page  is  scanned.  This  continues  until  either  a  favorable 
variable  is  found  or  all  basic  variables  have  been  priced.  The  list  of  basic 
variables  are  treated  in  a  wrap-around  fashion.  Let  £  denote  the  index  of  the 
most  recently  priced  entry,  with  £  initially  set  to  one.  The  fixed-page  pricing 
technique  can  be  described  mathematically  as  follows: 

Procedure  :  DJFixed-PageJPricing 
Input :  L 
Output :  £,  q. 
begin 

j  *  Q  * 

while(  j  <  |V|  &  q  =  0  )do 
i  <-  1;  p  <-  min{p,  |V|}; 
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while(  i  <  p  )d° 

k  «-  A{1)\ 

h-Xk 

Ak  *—  { 


/  *  price  p  basic  arcs  *  / 
j  *  get  basic  arc  index  associated  with  node  l  *  / 

t*  if  xk  <h\  . 

Xk-uk  if  Xk  >  Uk ;  /  *  get  the  bound  violation  value  *  / 

0  otherwise; 

max  r  A  A  -» 

{A*;,  Ag}; 

mod(^,  |V|)  +  1; 


argmax 

j  +  1 ;  ^ 


i  i  +  1;  j 

endwhile 

endwhile 

end. 

The  fixed-page  pricing  strategy  can  be  enhanced  by  maintaining  multiple  pric¬ 
ing  candidates.  A  candidate-queue,  which  is  a  cyclic  list  of  varying  ength  is 
ustd  to  keep  the  indices  of  the  favorable  variables  (  variables  violating  their 
bounds  ).  The  candidate-queue  pricing  technique  starts  with  a  general  scan 
pricing  a  set  of  basic  variables  to  select  the  most  favorable  variable  (  the 
variable  having  the  maximum  bound  violation  ),  while  placing  other  favora*> 
variables  into  the  queue.  At  each  iteration  a  subset  of  the  queue  entries  are 
priced  The  most  favorable  variable  becomes  the  leaving  variable,  and  other 
variables  that  price  favorable  are  placed  back  in  the  queue.  A  general  scan  is 
used  to  replenish  an  exhausted  queue,  where  each  scan  starts  where  the  previou 

one  has  ended. 

Let  Q,  called  the  candidate-queue,  be  a  sequence  of  entries  and  let t «  denote  the 
maximum  allowable  size  of  the  candidate-queue.  Let  Q[l]  denote  the  first  entry 
of  the  candidate-queue,  and  let  Q[i ... j ]  denote  entries  t  through  j  Let  Q&&W 
denote  the  concatenation  of  (<)  to  the  end  of  Q.  Let  p  denote  the  number  of 
queue  candidates  to  be  priced  in  each  iteration.  Let  l  denote  the 
node  associated  with  the  most  recently  priced  basic  variable  with  l  initialized 
to  one.  The  candidate-queue  dual  pricing  technique  can  be  represented  as 

follows: 

Procedure  :  D  -Candidate -Queue -Pricing 
Input  :£,Q. 

Output :  t, 
begin 


q  i -  +  0, 

while(  Q  ^  &  g  =  0  )do 

i  <-  1;  p  *-min{p,  |Q|}; 
while(  i  <  p  )do 

j*-Q[ l];  Q«-Q[2...|Q|]; 


/  *  price  p  queue  entries  *  / 
/  *  get  the  first  queue  entry  *  / 
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k  < A(  j)'  /  *  get  basic  arc  index  associated  with  node  j  *  / 

{h-x k  if  <  lk]  ,  .  .  ,  / 

xk-uk  if  xk>uk-,  /*  get  the  bound  violation  value*/ 

0  otherwise; 

f  f  I  _  “sm“  {An,  A,}; 

I  k  I  atgmm 

if  (  Afc  >  0  )then  Q  Q  &&  (k)\ 

i  <—  i  +  l; 

endwhile 

SJd,W=  o' )then  D.Generalurcantf, ,,Q);  /  »  repl»»h  the  «mply  queue  .  / 

The  mechanism  for  replenishing  an  empty  queue,  is  to  perform  a  general  pnang 
scanTn  *  basic  variables,  and  place  the  favorable  arcs  mto  the  queue.  The 
general  scan  procedure  can  be  stated  as  follows: 

Procedure:  D.Generalscan 
Input :  £,Q- 
Output :  £,  q,  Q- 
begin 

q  <—  Aq  <—  0;  j  <—  1; 
while(  j  <  |V|  &  q  =  0  )do 
i  < —  1; 

while(  i  <  k  &  |Q|  < w  )do  .  ,  j  p *  / 

k*-A(ty,  /  *  Set  basic  arc  mdex  associated  Wlth  nod  i  ‘ 

{lk  -  ik  if  Xk  <  h;  ...  ,  , 

Xk-Uk  \ixk>Uk\  /  *  get  the  bound  violatxon  value  *  / 

0  otherwise; 

r  g  1  argmax  {A  A  }; 

I  k  I  argmin  1  9 

iff  At  >  0  )then  Q  <—  Q  &&  (&); 

ih-iH- 1; J  *-j  +  mod(£,  |V|)  +  1; 

endwhile 

endwhile 

end. 

The  above  pricing  strategies  (maximum-infeasibility,  fixed-page  candidate- 

queueT 2  Lilafto  strategies  that  have  been  used  in  the  pnma 

For  problems  having  4000  or  more  nodes,  the  pricing  strategy  can  affect  the 

performance  of  the  dual  algorithm. 
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2.3  Steepest  Edge  Dual  Simplex  Algorithm 

In  the  simplex  method,  we  proceed  from  one  vertex  to  an  adjacent  one  along 
a  selected  edge.  When  using  standard  pricing,  the  edge  is  selected  such  that 
the  objective  function  improves  the  most  along  that  edge.  In  the  standard  dual 
simplex  method,  the  size  of  the  bound  violation  is  taken  directly  to  indicate 
the  gradient  improvement.  This  can  be  regarded  as  measuring  the  gradient  in 
the  framework  of  the  current  basic  variables.  In  the  steepest  edge  method,  the 
selected  edge  is  the  one  along  which  the  improvement  of  the  objective  function 
is  the  highest  with  respect  to  the  original  solution  space.  This  involves  taking 
into  account  the  variables  that  change  along  the  edge.  That  is,  the  selected 
edge  is  the  one  that  makes  the  largest,  most  obtuse  angle  with  the  gradient  of 
the  objective  function. 

It  is  easy  to  show  that  the  edge  directions  for  a  particular  basic  dual  solution 
are  given  by  the  rows  of  the  corresponding  basis  inverse,  B”1.  To  select  edges 
that  are  row-length  independent,  edge  normalization  was  introduced,  where 
weighting  factors  are  used  for  estimating  the  edge-length.  The  steepest-edge 
pricing  is  one  of  a  variety  of  methods  for  normalized  pricing,  in  which  the  bound 
violations  are  row  scaled.  However,  it  is  not  practical  to  calculate  the  edge- 
lengths  at  each  iteration.  Harris  [14],  used  a  dynamic  subset  of  the  original 
solution  space  as  a  working  framework.  An  approximation  of  the  weighting 
factors  are  obtained  using  the  Euclidean  norm  of  the  sub-vector  consisting  of 
just  the  components  in  the  current  working  framework.  That  is,  the  Harris5 
variant  of  the  steepest  edge,  called  Devex ,  chooses  the  edges  that  are  only  ap¬ 
proximately  steepest.  Goldfarb  and  Reid  [13],  developed  a  recurrence  formula 
for  the  squares  of  the  Euclidean  norms  of  the  edge  direction  when  pivoting 
from  one  feasible  basis  to  a  new  one.  Forrest  and  Goldfarb  [9],  introduced 
several  new  steepest  edge  algorithms  for  both  the  primal  and  the  dual  simplex 
methods. 


Let  rji  =  ||e,B*“1||2, i  =  l,...,m,  denote  the  length  squared  of  row  i  of  the 
basis  inverse.  Following  the  previously  described  dual  simplex  terminology,  the 
steepest  edge  dual  simplex  algorithm,  using  the  Goldfarb  and  Reid  recurrence 
formula,  for  a  given  dual  feasible  solution  defined  by  and  Afu  can  be 

stated  as  follows: 

Procedure  :  Dual  Steepest. Edge 
Input :  c,G,/,r,  u,B,A/l,A4* 

Output :  xB ,  xN ,  v* . 
begin 

7T  cbB~x ;  0-fc  <—  c*  -  7 vgk,  k  £ 
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r  <—  r  -  Y,k e/S,  9kh  ~  9kuk',  xB  B 

v* «-  EteM.  CkUk  +^keBckxk> 

for(  i  6  B  )do  z  <—  e‘B_1;  r?,-  *-  zz\  /  *  initiate  the  row  norms  *  / 

while(  3i  G  B  :  (it  <  h)  or  (i>  >  u,-)  )do 


argmax 


lA(k)-*lc 


if  xf  <  Ik 
if  if  >  uk 


|  fi^SM  if  if  >  Uk  j 

q  *-  A(g);  z  <-  e?‘5-1; 
w  *—  sign(u5  —  Xq)zG; 

£  _  {j  g  Xi  :  toy  <  0}  U  {j  G  Xu  :  u>y  >  0}; 
5  I  <_  naaxjgf  f  £i.  1  . 
s  J  argmaxJ  gf  l  w>  i  ’ 

y  «-  V;  V  =  {i:  yi  #  0}; 
f  (ZszM  if  *,</,  I 

LEsffA  if  Xq  >  Uq  J  ’ 

xf  <—  xf  -  A yiJ  eV]  xs  «—  x5  -f  A; 

f  .  ,,  f  <5  ifx^: 

crjb  «-  crfc  -  6w;jb,  &  G  JV;  crg  <-  <  . 


k  =  1, . . . ,  m; 


/  *  updated  row  *  / 

/  *  Ratio  Test  *  / 
/  *  updated  column  *  / 

/  *  flow  update  *  / 


<5  if 


— 5  if  =  uq  j  9 

v*  v*  +  cr5  A;  /  *  update  the  objective  value  *  / 

7T « —  7r  —  (5z;  t  <—  B_12; 

for(  i  =  1, . . . ,  m  :  i  #  g  )do  /  *  update  the  norms  *  / 

7/j  *-  max  (»Ji,  1  +  (^)  ^  ; 
endfor  ^ 

Vi  *-  (jfc)  w> 

B  <-BU  {s}\{qy,  tfu*-K\{»Y,  /*  basis  update  *  / 

if(  xq  =  lq  )then  X;  <—  M  U  {<?}; 
if(  i?  =  Uq  )then  Xu  <—  A4  U  {g}; 

endwhile 

end.  •  •  ..... 

In  the  above  dual  steepest  edge  procedure,  the  standard  maximum-infeasibihty 
strategy  is  used  for  selecting  the  leaving  variable.  The  partial  pricing  techniques 
were  also  extended  to  be  used  in  variants  of  the  steepest  edge  algorithm. 
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2.4  An  Initial  Dual  Feasible  Solution 

A  solution  x ,  defined  by  Mv,  is  said  to  be  dual  feasible  if  the  correspond¬ 
ing  dual  variables  (7t,A,/z)  £  X>.  Thus  far,  we  assumed  that  the  dual  simplex 
algorithms  are  initiated  with  a  dual  feasible  solution.  The  strategy  needed  to 
obtain  a  starting  dual  feasible  solution  was  not  evident.  A  starting  basis  B 
for  the  generalized  network  problem  is  obtained  by  adding  an  artificial  variable 
(£i,  i  =  1, . . . ,  m)  with  zero  bounds,  to  each  of  the  constraints  Gx  =  r.  There¬ 
fore,  the  set  of  dual  feasible  solutions  VA  =  £  3?m,  Xx  £  5Rn,  A^  £  3im,  \ix  £ 

S  :  GT7r-f  A1  -nx  =  c*,  tt  +  A*  - =  <£,  Ax  >  0,/ti*  >  0,  Ay  > 

0,/iy  >  0}.  In  order  to  assure  that  T>A  is  equivalent  to  D,  the  artificial  cost 
vector  c**  is  set  to  zero.  Therefore,  for  this  starting  basis,  7T  =  0,  the  arc  reduced 
costs  are  simply  the  arc  costs,  and  any  feasible  solution  of  *2),  say  (tt,  XT}fix)i 
has  a  corresponding  feasible  solution  (7r,  A®,  A^, where  A^  —  //  =  7r.  In 
the  following  sections,  we  introduce  techniques  for  obtaining  a  dual  feasible 
solution  from  any  starting  basis. 

Big-M  Method 

In  this  approach,  the  initial  dual  feasible  solution  is  obtained  by  setting  the 
original  variables  (zrt,i  =  l,...,n)  to  be  nonbasic  at  the  appropriate  bound 
based  on  the  sign  of  the  corresponding  reduced  costs.  That  is,  Mi  =  {j  £ 
B  :  (Tj  >  0},  and  Mu  =  {j  £  B  :  <Tj  <  0}.  Then,  the  artificial  variables 
£  =  f,  where  r  is  the  updated  right-hand-side.  An  infinite  upper  bound  is 
assigned  to  the  original  unbounded  variables,  where  infinity  is  approximated 
by  a  large  positive  constant,  usually  called  big-M.  The  dual  problem  becomes 
max(7r>x,/j)e'ZM  {rir  4  /A  —  ufj,}  ,  and  optimality  is  achieved  when  all  variables 
satisfy  primal  feasibility.  Since  the  artificial  variables  have  zero  bounds,  the 
corresponding  flows  will  be  forced  to  zero  at  optimality.  Also,  variables  with 
big-M  upper  bounds  will  be  driven  to  a  smaller  value.  Other  big-M  approaches 
have  been  suggested  to  obtain  a  dual  feasible  starting  basis  (  Murty  [20]  ). 

Two  Phase  Method 

The  big-M  approach  for  setting  the  initial  basis  is  not  recommended  because 
of  difficulties  in  making  the  appropriate  choice  of  M.  To  avoid  this  critical 
setting  of  the  big-M  value,  a  two  phase  approach  is  usually  used.  In  phase 
1  a  new  problem  associated  with  the  original  dual  problem  is  defined  whose 
optimal  solution  is  dual  feasible  for  the  original  problem.  The  initial  solution 
is  obtained  by  setting  the  original  variables  (a?,-,  i  =  1, . . . ,  n)  to  be  nonbasic  at 
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the  appropriate  bound  based  on  the  sign  of  the  corresponding  reduced  costs. 
Unbounded  variables  with  negative  reduced  costs  are  assigned  to  their  lower 
bounds,  resulting  in  a  dual  infeasibility  equal  to  the  sum  of  the  corresponding 
reduced  cost.  Therefore,  tfv  =  {j  £  B  :  a,  <  0,  uk  is  finite},  X  =  {j  g 
B  :  <Tj  <  0,u*  is  undefined},  and  A/}  =  {j  &  B  :  <Tj  >  0}  Ul.  The  phase 
1  objective  function  is  obtained  by  temporarily  setting  the  lower  and  upper 
bounds  to  zero  for  all  variables  and  introducing  a  new  right-hand-side,  fi  = 
Zk,j{-9k}-  Therefore,  the  artificial  variables  (  =  r,  and  the  phase  1  objective 
is  to  max^A^gtM  {fir}.  This  approach  is  similar  to  and  influenced  by  the  two 
phase  method  developed  by  Professor  Bixby  in  CPLEX  3.0.  Let  <Jj  denote 
the  current  reduced  cost  for  j  6  I.  The  contribution  to  dual  infeasibility  by 
xj  is  (0  -  (<Tj  -  6wj)).  Thus,  the  total  sum  of  infeasibility  is  J2jsi(5wi  ~ 
Cj).  Differentiating  with  respect  to  <5,  the  rate  of  change  of  the  sum  of  dual 

infeasibilities  is  ZjeiM  =  Eiei {-**0  =  *Eje rW‘>-  Let  *«  be  the 
leaving  variable.  Therefore,  the  rate  of  change  of  the  objective  function  is 
xq  =  e*x  =  e*B~lr  =  z(r  -  E ~  E SettinS  the  uPPer  and 
lower  bounds  on  all  variables  to  zero,  and  setting  the  right-hand-side  to  be 
defined  in  terms  of  the  variables  that  violate  dual  feasibility,  results  in  having 
the  rate  of  change  of  the  objective  function  to  be  the  same  as  the  rate  of 
change  of  the  sum  of  dual  infeasibilities.  This  approach  can  be  applied  to 
any  basis,  and  it  is  not  restricted  to  an  all  artificial  starting  basis.  Another 
intuitive  interpretation  of  this  setting  is  that  for  any  j  £  X  the  reduced  cost 
aj  =  Cj  -  ngj  <  0.  The  objective  function  to  drive  aj  to  be  nonnegative  is 
min  }.  Therefore,  for  all  the  dual  infeasible  variables  the  objective  function 
is  min{7rE jei93}- 

The  computational  steps  of  the  dual  simplex  algorithm  in  phase  1  are  modified 
to  incorporate  dual  infeasibilities.  The  pivoting  rules  are  as  follows: 


1.  the  leaving  variable  assumes  the  closest  of  the  upper  or  lower  bound,  and 
the  reduced  cost  will  assume  the  correct  sign.  Unbounded  variables  are 
allowed  to  move  only  to  the  lower  bound. 

2.  the  basic  variables  remain  dual  feasible  after  updating  the  basis. 

3.  the  entering  variable  assumes  a  dual  feasible  reduced  cost. 


In  this  approach  the  total  number  of  infeasibilities  is  reduced  without  restricting 
the  total  amount  of  infeasibility.  Other  methods  that  reduce  the  total  amount 
of  infeasibility  require  fewer  iterations,  but  need  some  extra  computational 
effort.  Belling-Seib  [4]  presents  similar  ideas  for  the  primal  simplex  method. 
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Following  the  previous  terminology,  the  phase  1  dual  simplex  algorithm  can  be 
stated  as  follows: 

Procedure  Dual-Phase  1 


Input  :  c,  G,  /,  r,  it,  M. 

Output :  x 3 ,  x v* . 
begin 

V  <Tk  *-  Ck  -  fgk,  k  6 

Afu  *—  {&  G  Af :  (Tk  <  0,  ujt  ^  oo};  X  * 
M  <—  {k  e  N  :  crk  >  0}  U  X; 


r  Eib6z{-5i}; 

while(  X  ^  <f>  )do 

q  <—  argmaxigV 


B~lr\ 


■  {fc  6  N  :  crjfe  <  0,  Wjt  =  oo}; 


f  — x®  if  x 

l  *£  ^ 5 


*f  <0 

B 


i2  f  ;  /  *  leaving  variable  *  / 

*  >  #  oo  »  '  &  ' 

g  4—  j4(£);  2:  4—  /  *  row  q  of  basis  inverse  *  / 

w  4—  sign(— x?)^G;  /  *  updated  row  *  / 

£  <-  {j  G  Aft  :  wj  <  0  ,  crj  >  0}  U  {i  G  Mi  :  uy  >  0  ,  a j  <  0}; 

if(  £  =  )then  £  <—  {j  €  I :  Wj  >  0}; 

maxje^  |  ZJ.  | ;  /  *  ratio  test  *  / 

/  *  updated  column  *  / 


5 
s 

V-B 


argmaxjEf 


if(  <  0  )then  M  «— M  U  {g}; 
if(  xq  >  0  )then  J\fu  4—  Nu  U  {$}; 

1  <- x?  -  Ayiy  i  €V; 


Vq 


A; 


/  *  flow  update  *  / 


if(  <5  >  0  )then  cr*  4—  <jfc  —  k  €  N;  <rq  sign(— £?)<5; 


B  4— jBU  {s}\{g};  M  4—  M\{s};  Nu  Mx\{s};  /*  basis  update  *  / 
for(  I:  e  I  &  <Tfc  >  0  )do 

«—  —  I?”1#*;  /  *  adjust  the  flow  *  / 

1 4—  J\{fc};  /  *  variable  fc  became  feasible  *  / 

endfor 
endwhile 


dual(c ,  G,  /,  r,  u,ByN,  xBy  xN ,  t?*); 

end. 

The  search  for  the  leaving  variable  in  phase  1  of  the  dual  simplex  method  is 
enhanced  by  exploiting  the  infrequency  (small  number)  of  the  dual  infeasibilities 
and  the  basis  structure  of  generalized  networks.  The  search  is  limited  to  small 
parts  of  the  basis  components  associated  with  the  nodes  in  X  =  {i  G  V  :  g\  ^ 
0 yj  G  X};  precisely  the  path  from  a  node  in  X  to  the  root  of  the  corresponding 
rooted- tree  component,  and  the  nodes  in  the  cycle  along  with  the  nodes  in  the 
path  from  a  node  in  X  to  the  cycle  of  the  corresponding  one-tree  component. 
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3  BENCHMARK  PROBLEMS  AND 
TESTING  ENVIRONMENT 

A  set  of  generalized  network  test  problems  were  randomly  generated  using  a 
modified  version  of  NETGEN  by  Klingman  et  al  [16].  The  modified  code,  gen¬ 
erates  a  variety  of  generalized  network  problems  based  on  the  user  controlled 
problem  description.  The  main  inputs  are  the  number  of  nodes,  the  number 
of  arcs,  the  number  of  source  nodes,  the  number  of  sink  nodes,  the  number  of 
transshipment  sources,  and  the  number  of  transshipment  sinks.  In  addition, 
the  minimum  supply  and  the  maximum  supply  specifications  bounds  the  total 
supply.  The  %  arcs  with  high  cost  parameter  specifies  the  number  of  arcs  hav¬ 
ing  the  maximum  cost,  and  the  remaining  arc  costs  are  randomly  generated  on 
a  user  specified  interval.  Similarly,  the  %  arcs  with  bounds  parameter  specifies 
the  percentage  of  bounded  arcs,  where  the  generated  arc  bounds  are  uniformly 
distributed  on  a  given  interval.  A  seed  is  supplied  to  a  random  number  gener¬ 
ator  that  allows  one  to  reproduce  the  problems. 

Our  modified  code  can  be  divided  into  two  main  segments.  In  the  first  segment, 
a  skeleton  network  that  guarantees  the  problem  connectivity  and  feasibility  is 
generated.  Initially,  the  total  supply  is  randomly  distributed  among  the  source 
nodes,  including  the  transshipment  source  nodes.  Then,  for  each  source  node 
a  chain  is  created  to  carry  the  flow  through  a  series  of  distinct  transshipment 
nodes,  to  be  randomly  allocated  among  a  randomly  chosen  subset  of  sink  nodes. 
These  chains  are  mutually  exclusive  except  for  the  sink  nodes.  In  the  second 
segment,  randomly  generated  arcs  are  appended  to  the  network  to  fill  out  the 
total  required  arcs.  The  final  network  will  contain  approximately  the  required 
number  of  arcs.  For  more  detail  about  generating  the  arc  costs  and  bounds  see 
Klingman  et  al  [16].  For  the  skeleton  arcs,  the  from  and  the  to  multipliers 
are  set  to  1  and  -1,  respectively.  For  each  of  the  remaining  arcs,  the  two 
multipliers  are  randomly  selected  from  the  continuous  interval  [—1,1),  while 
insuring  that  the  first  multiplier  is  nonzero. 

The  comparison  and  examination  of  the  efficiency  for  different  solution  methods 
is  of  practical  importance,  and  the  selection  of  the  test  set  used  in  the  com¬ 
parison  is  very  crucial.  Therefore,  we  generated  ten  test  problems  of  various 
sizes  and  structures.  The  number  of  nodes  ranged  from  1,000  to  15,000,  and 
the  required  number  of  arcs  ranged  from  3,000  to  90,000.  The  characteristics 
of  the  test  problems  are  given  in  Tables  1  and  2.  The  naming  convention  is 
GJhJa  for  a  generalized  problem  having  mK  nodes,  and  hK  arcs.  All  testing 
was  performed  on  a  60  MHz  DECStation  5000/260  running  the  Ultrix  4.3a  op¬ 
erating  system.  The  reported  performance  measures  are  the  number  of  pivots 
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Table  1  Characteristics  of  test  problems  one  through  five. 


Specification 

Problem 

G-01-03 

G-01-06 

G-02.06 

G-02.12 

G_04_12 

Seed 

85 

173 

319 

51 

920 

Nodes 

IK 

IK 

2K 

2K 

4K 

3K 

6K 

6K 

12K 

12K 

Sources 

100 

100 

200 

■ESI 

HI 

Sinks 

100 

100 

200 

mmm 

Trans.  Sources 

50 

■HI 

100 

100 

100 

Trans.  Sinks 

HHE39I 

100 

mem 

Chain  Length 

8 

io 

20 

25 

Min  Cost 

-IK 

-2K 

-IK 

-3K 

-4K 

Max  Cost 

10K 

10K 

10K 

10K 

HUI 

Min  Supply 

iHoa 

10K 

50K 

50K 

w*mm\ 

Max  Supply 

20K 

20K 

50K 

50K 

70K 

Min  Arc  Bounds 

IHHKEI 

10 

10 

mmm 

Max  Arc  Bounds 

IK 

IK 

IK 

IK 

IK 

High  Cost  Arcs 

58% 

40% 

29% 

■Hil 

Bounded  Arcs 

70% 

99% 

87% 

73% 

92% 
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Table  2  Characteristics  of  test  problems  six  through  ten. 


Specification 

Problem 

G-04.20 

G-09.20 

G-09.45 

G.15.45 

G-15-90 

Seed 

11 

194 

239 

357 

421 

Nodes 

4K 

9K 

wmm 

Arcs 

20K 

20K 

45K 

45K 

—sail 

Sources 

200 

200 

mmm 

^uuu 

Sinks 

■Mil 

200 

Ha 

W/ESMk 

2000 

Trans.  Sources 

100 

100 

250 

1000 

100 

100 

250 

1000 

251 

50 

■KjUgl 

-2K 

-llT 

mil 

10K 

30K 

H§3 

80K 

HU 

wmm\ 

■ 

90K 

BEsai 

1 

10 

100 

|  100 

IEBE3Eg5Ea 

1 

■Ha 

IH^ 

I9MKL9 

iheii 

jsnignsmiEtiM 

IHEI&J 

IHi 

!SH! 

SHHaffiMl 

1 

«Ksa 

mtw 

IMEzJ 

IHH1I 
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(  or  iterations  )  and  the  user  central  processing  time  rounded  to  two  deci¬ 
mal  places.  The  user  cpu  time  is  obtained  using  the  standard  Ultnx  function 
getrusage,  that  returns  information  describing  the  system  resource  utilization. 
This  information  includes  the  cumulative  user  cpu  time,  measured  in  microsec¬ 
onds,  consumed  to  the  point  of  the  function  call.  In  the  implementation  of  the 
simplex  algorithm,  getrusage  is  called  immediately  preceding  and  immediately 
following  the  call  to  the  optimizer.  The  user  cpu  time  consumed  is  obtained  by 
subtracting  the  first  call  returned  time  from  the  second  call  returned  time. 


4  COMPUTATIONAL  EXPERIENCE 
WITH  CPLEX  3.0 

To  provide  a  comparison  base,  the  test  problems  were  solved  using  the  four  op¬ 
timization  algorithms  available  in  CPLEX  3.0  (  a  state-of-the-art  solver  [24]. 
This  version  includes  efficient  implementations  of  the  primal  simplex,  the  dual 
simplex,  the  network  simplex,  and  the  barrier  algorithms.  All  CPLEX  runs 
were  made  using  the  default  parameter  settings,  which  activate  the  presolver 
and  the  aggregator  problem  preprocessors.  The  presolver  helps  simplify,  re¬ 
duce,  and  eliminate  redundancies  in  the  problem  presentation.  The  aggregator 
attempts  to  make  simplifying  substitutions  that  reduce  the  basis  size. 

In  the  CPLEX  dual  simplex  optimizer,  the  dual  problem  is  solved  using  the 
primal  simplex  presentation,  where  all  the  linear  algebra  is  carried  out  on  the 
associated  primal  basis.  The  network  optimizer,  recognizes  and  exploits  the 
embedded  network  structure  of  the  problem,  and  uses  an  efficient  specialized 
algorithm  on  the  pure  network  portion  to  obtain  an  advanced  start  for  the  dual 
optimizer.  The  barrier  optimizer,  is  an  efficient  implementation  of  the  primal- 
dual  method  fully  integrated  with  the  predictor  corrector  basis  crossover  of 

Megiddo  [17]. 

The  empirical  results  for  solving  the  test  problems  using  the  four  CPLEX 
solvers  are  presented  in  Table  3,  where  the  numbers  in  parenthesis  represent 
the  number  of  phase  I  pivots.  From  the  first  four  problems,  it  is  clear  that  the 
barrier  solver  performs  poorly  for  generalized  network  problems,  and  it  was  not 
run  for  the  remaining  ones.  The  CPLEX  dual  simplex  solver  is  approximate  y 
two  times  faster  than  the  primal  simplex  solver  for  all  problems  except  G.01.03, 
G  15  45  and  G-15J90.  The  over  all  time  performance  is  1.4  times  faster  than 
the  primal  algorithm.  Using  the  default  settings,  the  CPLEX  network  solver 
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failed  to  extract  any  network  component  for  three  of  the  test  problems,  and  it 
terminated  with  an  unbounded  network  for  problems  G-09.45  and  G-15-45. 


5  EMPIRICAL  ANALYSIS 

The  specialized  big-M  dual  simplex  algorithm  is  implemented  in  C,  and  the  per¬ 
formance  in  solving  the  ten  test  problems,  using  the  three  pricing  techniques, 
is  shown  in  Tables  4,  5,  and  6.  For  all  testing  M  is  set  to  109.  Comparing 
the  big-M  standard  (  non-normalized  )  dual  RAMSES,  using  the  three  pricing 
strategies,  and  the  CPLEX  3.0  dual  solver,  it  is  clear  that  the  dual  RAMSES 
with  a  candidate- queue  pricing  is  the  superior.  For  the  G-15-90,  it  performed 
approximately  30  times  faster  than  CPLEX,  13  times  faster  than  RAMSES 
with  a  maximum-infeasibility  pricing  strategy,  and  28%  faster  than  RAMSES 
with  a  fixed-page  pricing  strategy.  On  average,  for  the  ten  problems  it  per¬ 
formed  20  times  faster  than  CPLEX,  9  times  faster  than  RAMSES  with  the 
maximum-infeasibility  pricing  strategy,  and  25%  faster  than  RAMSES  with 
fixed-page  pricing  strategy.  The  fixed-page  pricing  strategy  resulted  in  only 
a  2%  increase  in  the  total  number  of  pivots,  while  the  candidate-queue  pric¬ 
ing  strategy  resulted  in  a  19%  increase  in  the  total  number  of  pivots  over  the 
maximum-infeasibility  pricing  strategy. 

For  the  three  pricing  strategies,  the  standard  (  non-normalized )  RAMSES  per¬ 
formed  consistently  much  faster,  and  required  fewer  pivots  than  the  steepest 
edge  RAMSES.  This  can  be  ascribed  to  the  generated  problem  characteris¬ 
tics  of  having  arc  multipliers  in  the  range  [-1,1],  and  uniformly  distributing 
the  arcs  among  the  network  nodes,  which  results  in  basis  inverse  rows  having 
approximately  the  same  length. 

The  specialized  two  phase  dual  simplex  algorithm  is  implemented  in  C  and  the 
performance  in  solving  the  ten  test  problems,  using  the  three  pricing  techniques, 
is  shown  in  Tables  7,  8,  and  9.  It  is  clear  that  the  two  phase  dual  RAMSES 
with  a  candidate-queue  pricing  strategy  is  the  superior  implementation.  For 
the  G-15-90,  it  performed  approximately  31  times  faster  than  CPLEX,  13  times 
faster  than  RAMSES  with  a  maximum-infeasibility  pricing  strategy,  and  29% 
faster  than  RAMSES  with  a  fixed-page  pricing  strategy.  On  average,  for  all 
test  problems  it  performed  22  times  faster  than  CPLEX,  10  times  faster  than 
RAMSES  with  a  maximum-infeasibility  pricing  strategy,  and  22%  faster  than 
RAMSES  with  a  fixed-page  pricing  strategy. 


Table  3  Solving  generalized  networks  using  CPLEX  version  3.0 
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Table  6  Results  for  big-M  dual  RAMSES  with  candidate-queue  pricing. 
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Table  7  Results  for  two  phase  dual  RAMSES  with  maximum-infeasibility  pricing. 
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Table  8  Results  for  two  phase  dual  RAMSES  with  fixed-page  pricing. 
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Chapter  6 


6  SUMMARY  AND  CONCLUSIONS 

Generalized  networks  are  used  to  model  a  diversity  of  practical  problems  that 
include  class  scheduling,  machine  loading,  manpower  planning,  resource  distri¬ 
bution  systems,  and  currency  exchange.  This  paper  describes  the  development 
of  efficient  techniques  for  solving  generalized  network  problems.  As  in  the  pri¬ 
mal  simplex  algorithm,  the  dual  simplex  method  is  initiated  by  an  all  artificial 
basis.  In  Section  2.4  we  introduce  two  approaches  to  obtain  a  dual  feasible  solu¬ 
tion  from  any  starting  basis.  In  the  two  approaches,  the  nonbasic  variables  are 
assigned  to  the  appropriate  bound  based  on  the  the  sign  of  the  corresponding 
reduced  costs.  In  the  first  approach,  a  big-M  value  is  used  to  approximate  the 
upper  bound  for  unbounded  variables.  For  the  three  dual  pricing  strategies, 
the  big-M  dual  RAMSES  code  with  a  candidate-queue  is  superior.  It  achieves 
a  performance  30  times  faster  than  the  CPLEX  3.0  dual  optimizer  for  problem 
G-15-90,  and  20  times  faster  for  all  the  benchmark  problems.  However,  the 
success  of  the  big-M  method  is  contingent  upon  the  selection  of  a  proper  big-M 
value.  In  the  second  approach,  this  critical  selection  is  avoided  by  using  a  two 
phase  method,  where  dual  infeasibilities  are  temporarily  permitted  in  phase  1 
for  unbounded  variables  with  negative  reduced  cost.  For  the  three  dual  pricing 
strategies,  the  two  phase  dual  algorithm  with  a  candidate-queue  is  superior. 
It  achieves  a  performance  31  times  faster  than  the  CPLEX  dual  optimizer  for 
problem  GJ.5-90,  and  22  times  faster  for  the  ten  test  problems. 
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